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SUMBfARY 


This report describes the formulation and development of a computer 
analysis for the calculation of streamlines and pressure distributions around 
two-dimensional (planar and axisymmetric) isolated nacelles at transonic 
speeds. The computerized flow field analysis is designed to predict the 
transonic flow around long and short high-bypass-ratio fan duct nacelles with 
inlet flows and with exhaust flows having appropriate aerothermodynamic pro- 
perties. The flow field boundaries are located as far upstream and down- 
stream as necessary to obtain minimum disturbances at the boundary. The far- 
field lateral flow field boundary is analytically defined to exactly represent 
free-f light conditions or solid wind tunnel wall effects. 

The inviscid solution technique is based on a Streamtube Curvature 
Analysis. The computer program utilizes an automatic grid refinement proce- 
dure and solves the flow field equations with a matrix relaxation technique. 

The boundary layer displacement effects and the onset of turbulent separation 
are included, based on the compressible turbulent boundary layer solution 
method of Stratford and Beavers and on the turbulent separation prediction 
method of Stratford. 

This computer program has the capability of calculating the pressure dis- 
tributions and flow fields, Including viscous displacement effects, on a variety 
of internal and external shapes. The location of incipient turbulent boundary 
layer separation is identified, if and when the calculated pressure gradients 
are sufficient to cause it. The computing times are relatively short (2-6 
minutes on a CDC 6600) depending on the complexity of the problem. The pre- 
dicted pressure distributions have been compared with the through-flow nacelle 
test results from the NASA-Langley 16- foot tunnel. 
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ANALYTICAL METHOD FOR PREDICTING THE PRESSURE DISTRIBUTION ABOUT A NACELLE 

AT TRANSONIC SPEEDS 

PART I - STREAMTUBE CURVATURE ANALYSIS 
J.S, Keith, D.R. Ferguson, C,L. Merkle, and P.H, Heck 
General Electric, Aircraft Engine Group 
Evendale, Ohio 45215 


1.0 INTRODUCTION 


Aircraft are being designed with the NASA-developed supercritical wing 
to fly at cruise Mach numbers approaching one. The need for low- installed 
drag and high-drag-divergence Mach number nacelle installations is extremely 
critical to the success of this design. Design techniques are required to 
evaluate these nacelles on an isolated basis and then on an installed or 
integrated basis. 

For this reason NASA has begun a program to provide design information 
for low-drag, high-drag-divergence Mach number isolated nacelles suitable 
for use with advanced high-bypass- rat io turbofan engines. One element of 
such a program is the development of a method to predict the inviscid pres- 
sure distribution and flow field about an arbitrary axisymmetric ducted body 
at transonic speeds. The prediction technique will provide the means to 
conduct parametric studies so that the nacelle design criteria could be 
evaluated to select configurations for further experimental investigations. 
The prediction technique would provide guidance during wind tunnel testing 
to develop nacelle shapes which would minimize drag within given design 
restraints. 

Several techniques of solving the inviscid equations of motion about 
arbitrary two- or three-dimensional bodies at transonic speeds are presently 
available; however, there are no ccxnputer programs available which treat air 
inlet or nacelle configurations. The objective of the development of this 
computer analysis was the prediction of flow fields about isolated nacelles 
at transonic conditions. The solution technique was further specified to 
give accurate results consistent with the requirement of relatively short 
computing time per input case as compared to that required for a time- 
dependent finite difference method of solution. The method utilized to com- 
pute the flow field is the Streamtube Curvature Relaxation technique. 

The Streamtube Curvature Method (STC) of solving planar and axisym- 
metric external flows has not been discussed significantly in the literature; 
however, the method is a very natural one. For example, engineers frequently 
rely on one-dimensional compressible flow relationships for a first-order 
solution to ducted flows. The STC approach is similar except that a number 
of confluent streamtubes, with slightly different properties, are added to- 
gether to obtain the total flow in the channel. Each streamtube is handled 
in much the same way as is the one streamtube in the one-dimensional problem. 
In the limit, as the size of the individual streamtubes approaches zero, the 
STC method satisfies the inviscid equations of motion exactly. 


This report describes the method of analysis used to apply the Stream- 
tube Curvature Relaxation technique, the numerical procedure for computer- 
ization of the analysis, and examples of correlations of predicted flow 
fields on nacelles at transonic speeds with wind tunnel test data. 

The computer program source deck, together with a user’s manual, is 
available from COSMIC (Computer Softwear and Information Center), Burrows 
Hall, University of Georgia, Athens, Georgia 30601. The program is written 
in CDC Fortran 2.3 source language, except for three subroutines in Compose 
1.1 language. The computer program has been checked out for the CDC 6600 
machine. 
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2.0 SELECTION OF THE METHOD 


Known methods for solving transonic flow fields may be divided into two 
categories — time- dependent and iterative. Time-dependent methods have 
achieved much popularity because both the subsonic and supersonic portions 
of the flow field, in most cases, are solved by the same algorithm. Thus, 
with a rather simple calculating procedure a difficult mathematical problem 
is computed. In the iterative method, however, the calculation formula must 
reflect the mathematical nature of the equation and a switching, depending 
upon Mach number, to the appropriate formula is required at each calculation 
point. It is in this way that the different physical characteristics of the 
subsonic and supersonic regions come into play. 

Iterative methods are quite new. To the authors' knowledge, the first 
demonstration of a general, numerically consistent, iterative method for 
solving transonic flows occurred in 1970 (ref. 1). This was the small per- 
turbation method of Murman and Cole for flow past airfoils without lift. 
Recently, extensions to the method have been presented by Stegger and Lomax 
(ref. 2). 

Although there undoubtedly are many variations, we may think of an 
Iterative method as one in which the equation for the unknown fluid d 3 mamic 
property at each of the net points is solved by: (1) writing a linear ap- 
proximation to this equation and, (2) solving the resulting system of equa- 
tions simultaneously. Because of the linear approximation, this process is 
repeated several times (say 3 to 10) before convergence is obtained. 

In contrast to solving the field simultaneously, time- dependent methods 
compute the wave motion of a disturbance as it travels from one part of the 
flow field to the other. A steady- state result is obtained only after all 
wave reflections have dissipated to a relatively small level. Although the 
time-dependent method of updating the flow properties can be likened to an 
iteration process, clearly the most rapid solution will be obtained when 
the flow field variables are all corrected simultaneously and when this cor- 
rection is not limited by (computationally) slow wave transits. Therefore, 
as a rapid analysis tool, the iterative method is most attractive. 

Of the many different representations of the fluid dynamic equations, 
the number which can be solved by the iterative method across the transonic 
region are, perhaps, limited. Here the simplest and most general forms of 
the equations are chosen, namely, those which apply along streamlines (f = 
constant lines) and those which apply along lines which are orthogonal to 
the streamlines (§ s constant lines). 

Across the streamlines, the continuity and momentum equations are: 



(§ = Const) (1) 
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Momentum: 


(a) Normal form: 


t 


(b) Crocco form: 


(| = Const) (2a) 


1 a^=.cv2*|5.Tf 

2 dn an an 


(§ = Const) 
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Along the streamlines the energy and momentum equations are: 
Momentum : 

DS 


Energy ; 


Ds 


m 

Ds 


= 0 


= 0 


(Y = Const) (3) 


(Y = Const) (4) 


where, the independent variables s and n are the distances measured along 
and across the streamlines, respectively. 

The solution method is an extension of the streamline curvature method. 
It may be briefly described as follows: First, a crude grid of streamlines 

and orthogonal lines is assumed (refer to Figure 1); second, the curvature 
of the streamlines at each of the grid points is evaluated; third, the mo-^ 
mentum equation is integrated along a line normal to the streamlines to 
obtain velocity, and the continuity equation is integrated to determine the 
’’correct’* streamline positions (for the assumed curvature field). These 
are indicated by the ”x** in Figure 1. Fourth, an adjustment (6n) is com- 
puted by considering: (1) the difference between the computed and assumed 

streamline positions and, (2) the effect of the implied curvature modifica- 
tion in the integrated momentum equation. Finally, the streamlines are re- 
positioned by the 6n values. 

Because the movement of any one grid point alters the velocity at 
nearby points through a change in curvature, it is highly desirable to 
account for these interrelating point adjustments simultaneously. The 
utilization of a simultaneous solution procedure, employed here, is not 
part of the classical streamline curvature method (refs. 3, 4, 5), In 
comparison, the classical method yields calculation times which are very 
slow, especially for a closely spaced calculation grid. In concept, the set 
of simultaneous equations for the normal streamline adjustments is formulated 
from the finite difference equivalent of the following equation: 
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1. Assume a Crude Grid 



2. Evaluate Curvature 


3. Integrate the Cross-Stream 
Momentum Equation and the 
Continuity Equation to Determine 
the "Correct" Streamline Positions. 



4. Solve the Matrix Equation for 
6n and Move the Grid Points. 


Figure 1. Solution Technique. 


S 


( 5 ) 


5^(6n) (l-M^) 9^(6n) _ „ 

where : 

6n = Required streamline adjustment in the normal direction 
Y = Stream function 

s = Curvilinear distance along a given streamline 

M = Mach number 

pV = Flow per unit area 

F = Driving (or error) function derived from the solution to 
the integral continuity and normal momentum equations 

This equation is derived in Appendix A for the special case of isen- 
tropic two-dimensional flow. (These limiting assumptions are utilized only 
to maintain simplicity of illustration; they are not part of the computer 
program.) From a mathematical point of view, the above equation is similar 
to the small perturbation form of the velocity potential equation employed 
by Murman and Cole (ref. 1). 

+ a-„2, l2| = 0 (6) 


V 

a 


« 1 , 



In either case, it is possible to numerically solve the equations for either 
subsonic flow or supersonic flow by changing the finite difference star from 
a subsonic representation to a supersonic representation as illustrated in 
Figure 2, Notice that the supersonic representation includes no points down- 
stream of the cross-stream line, reflecting the physical reality that dis- 
turbances downstream will not be felt upstream. The star- switching process 
is directly related to the coefficient (l-M^); and, because this coefficient 
is zero at unity Mach number, the switch from one star to the other is per- 
formed smoothly. 

The extended streamline curvature method, here referred to as the 
Streamtube Curvature (STC) method, appears to have the advantage that it is 
applicable to the calculation of nonsmall perturbation transonic flows. 
Considerable complexity is introduced when Equation 6 is expanded to allow 
the vertical component of velocity, v, to be the order of magnitude of the 
axial component u. In this case when the grid system is not aligned with 
the flow direction, a cross-derivative term: 
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Figure 2. 


Finite Difference Stars for Subsonic 
and Supersonic Flow. 



uv 5^^ 

Sx9y 

appears in the differential equation, and the star-switching concept (as 
explained above) cannot be applied. On the other hand, with the intrinsic 
coordinate system utilized in the STC procedure, the arms of the star are 
always oriented in the streamwise and cross- stream directions, and the star- 
switching algorithm is always appropriate. 

Star switching is one of the requirements for the numerical solution 
to either Equations 5 or 6 in a mixed subsonic and supersonic region. A 
second requirement is a reasonably accurate evaluation of the nonlinear 
coefficient term (1-M^). In this respect, the velocity potential method 
is superior. During the current contract, it has been found that the coef- 
ficient term in the differential equation (Equation 5 above) cannot readily 
be evaluated to the accuracy necessary for a convergent transonic solution 
when the streamwise spacing between the calculation points is too small. 

This evidently is due to the method used for evaluating the velocity (and 
Mach number) at field mesh points. The velocity calculation, using Equa- 
tion 2b, requires the evaluation of curvature by a second- order numerical 
differentiation which is subject to large errors (order of 1/A3^) when the 
streamwise spacing (AS) is small. 

Aside from the limitation just cited, the streamtube curvature method 
is extremely powerful as indicated by the following features; 

• No additional complexities arise when the flow is rotational. 

• The slip line between an exhaust jet and the external flow can 

be handled precisely. (The procedure is to consider two coincident 
streamlines. Their position and pressure are the same; their 
velocity and stagnation properties may be different,) 

• The streamline/orthogonal line grid provides a mapping of the flow 
field into a rectangular domain. This is helpful from the stand- 
point of ccHTiputer program organization. 

The STC Program has also been designed to: 

• Handle multiple streams 

• Place grid points at locations in the flow field where they are 
needed, as determined by local variations of the dependent 
variables 

• Allow external flow analysis by incorporating matched near-field 
and fai^field solutions. The far-field solutions are obtained 
analytically utilizing small perturbation theory. 


8 



3 . 0 AN OUTLINE OF THE CALCULATION STEPS 


The operations performed by the STC Program may be outlined as follows: 

1. Define the flow regions and locate (approximately) the "primary" 
orthogonals and the streamlines which divide the internal and 
external flows, 

2. Refine the grid as required by inserting additional streamlines 
and orthogonal lines between those already existing. 

3. Compute the streamline angles and curvatures. 

4. Compute the orthogonal line angles and move the grid points along 
the streamlines to obtain orthogonality, 

5. Compute the velocities on the "far-field" boundary, 

6. Adjust the flow rates in the exhaust streams, if any, to meet the 
calculated choking flow rate. 

7. Integrate along each orthogonal the mcmientum and continuity 
equations (Equations 1 and 2). 

8. Detemlne if the streamline positions are within a "rough tolerance," 
If so, return to Step 2 for additional grid refinement (unless grid 
refinement limits have already been reached). Otherwise, continue 

to Step 9. 

9. Determine if the streamline positions are within final tolerance. 

If so. Jump to Step 13. Otherwise continue to Step 10. 

10. Set up the matrix equation for the streamline correction, 6n. 

11, Solve the matrix equation. 

12. Modify the streamline positions by 6n, and return to Step 3. 

13, Calculate and print the output quantities; then return to Step 1 
for the next case, if any. 

The first operation includes reading the card input lor a description 
of the gecmietry and flow properties. The computer program has been written 
to have general capability lor analyzing a great variety of configurations. 

The first step in the programmed logic is to develop a table of orthogonals 
or calculation stations for each of the several flow regions. The regions 
are determined as illustrated in Figure 3, so the calculation can proceed 
from upstream to downstream. The boundary of each region is defined as a 
primary orthogonal. As shown in Figure 4, the initial grid which is de- 
veloped contains only the primary orthogonals and the double streamlines 
which separate the various streams. 


9 








The second step in the computational procedure is the grid refinement. 
The very crude grid, obtained in Step 1, is refined before the first solution 
of the flow field equations is executed, A new orthogonal is placed within 
each region and, likewise, a streamline is inserted in the middle of each 
channel. In the external channel, additional streamlines are placed close 
to the body. After the solution has been obtained for this net, the grid 
intervals are halved as required. This may be likened to the steps taken 
when one **flux plots” a flow field by hand. First, major flow lines and 
nomals are sketched in, and then more and more streamlines and orthogonal 
lines are added until the desired resolution is obtained. At each step in 
the process, the positions of the lines are adjusted to meet the correct 
solution requirements. The procedure automatically provides for grid re- 
finement in regions of high curvature and high acceleration or deceleration. 
The streamline and orthogonal lines which are added between existing lines 
are not required to span the field if only local refinement, near the body, 
is required. The refinement procedure presently built into the program uses 
a criteria involving the distance and velocity increment between grid points. 
These refinement criteria are discussed in detail in Reference 16, 

The third step in the method is to determine the angles and curvatures 
of the streamlines at each grid point. For subsonic portions of the flow 
field, this is performed by fitting a piecewise continuous cubic polynominal 
in a coordinate system which is locally rotated for each interval. The re- 
sulting fit is analogous to the curve produced by a beam which is loaded by 
discrete forces to pass through the given grid points. The locally rotated 
coordinate system removes the restriction that requires the slope to be 
small. For grid points located in a supersonic region, backward difference 
formulas are employed. Either 3- point or 4- point formulas may be 
optionally selected. Again the coordinate system is rotated so that slopes 
in the curve- fitting coordinate system are small. 

In the fourth step, the orthogonality of the grid points is checked and 
points are moved along the spline curve as required to achieve normal in- 
tersections between the two sets of lines. Also, the normal distance, n, 
is computed for each grid point as measured from the lower boundary of the 
orthogonal. 

When the initial grid is set up, a boundary is placed some distance 
away from the body. This boundary becomes the interface between the near- 
field and the far- field solutions. The near- field is computed by the stream- 
tube curvature method, and the far-field is computed by linear small per- 
turbation theory. In the process of iterating, this boundary (which is 
also a streamline) will float so that its shape and velocity distribution 
are matched by both the inner and outer solutions. In practice, the shape 
of the interface streamline (also referred to as the far-field boundary) is 
first assumed. Using the far-field equations, the velocity distribution is 
calculated. This is Step 5. These velocities are subsequently employed in 
the near-field analysis and from this comes a revised shape for the inter- 
face streamline. Revised velocities will then be computed in Step 5 during 
the following iteration cycle, and so forth. 
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step 6 is the modification, as required, of the flow rates of the ex- 
haust streams. For boattail analysis of nacelles, the Internal geometry 
of the exhaust passage is required input to the STC program. Because of 
streamline curvature effects, the discharge coefficient for the nozzle will 
be somewhat less than unity. The user, however, may input a flow rate based 
on unity discharge flow coefficient or, for that matter, any approximate 
value. Determination of the velocity distribution across the throat of the 
nozzle will be determined within the STC framework, and the evaluation of 
the maximum "choked" flow rate is Step 6 of the calculation procedure. 


Step 7 is the solution of the flow field equations per se. This section 
of the program is referred to as the "flow balance;" Equations 1 and 2 are 
Integrated. In the external regions of the field, the momentum equation is 
Integrated from the far- field Interface boundary to the body (or to the 
centerline or lower boundary, whichever exists). The integral form of the 
momentum equation is: 


InVfcZ . inv^ 


2 I Cdn 


(7) 


where: 

~ Velocity as determined in Step 5 along the far-fleld streamline 

s Velocity at any streamline with orthogonal distance nj^ 

n^g = Distance measured along the orthogonal to the far-fleld 
streamline 


Although not reflected in Equation 7, the effect of varying total pres- 
sure behind a shock wave is also Included. The method for handling this is 
presented in Section 5.5. Also, if a slip-line occurs in the field, the ve- 
locity Jump equations 



are employed where the subscripts (+) and (-) denote conditions on the 
streamlines above and below the slip-line, respectively. 
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The velocity, total temperature, and total pressure allow determination 
of the density at each grid point , and the inverse product of density and 
velocity is integrated to find flow area. 

|v 

The cumulative flow areas calculated by Equation 9 are compared with 
the geometric areas of the streamlines used in Step 3. The difference be- 
tween these two values is used as a convergence check (Steps 8 and 9) and 
in the streamline correction equation. Step 10, 

For internal flow orthogonals, the velocity at the outer boundary 
(Vff in Equation 7) is not known. Instead, an iteration process is employed 
whereby the outer boundary velocity is varied to obtain a match in the cal- 
culated geometric passage area. 

In Steps 10 and 11, the proper adjustment of the streamline positions 
is determined; and, in Step 12, the grid points are moved in the normal 
direction by this computed adjustment. 

The iterative sequence is to start with a crude grid, as noted above, 
and to repeat Steps 3 through 12 until the flow balance error is small. 

This is often accomplished in one or two iterations. The grid is then re- 
fined to the next level, and the field is reconverged. The ref inement /con- 
vergence process is continued until the grid refinement criteria is satis- 
fied, or alternately, until computer storage limits are reached. At this 
point, additional loops through Steps 3 to 12 may be performed until the 
flow balance error is satisfactory. 

In the next section, examples are shown of STC predictions; and, in 
Section 5, the details of the numerical pi^ocedures are presented. 
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4.0 EXAMPLE RESULTS 


4.1 SHORT DUCT FAN INSTALLATION 

The development of this computer analysis has been addressed to tran- 
sonic pressure distributions on typical Jet engine installations in isolated 
nacelles. The analysis was required to handle inlet flows and exhaust flows 
with correct aerothermodynamic properties. An example which demonstrates 
the capabilities is the short duct fan installation. 

The installation and the predicted flow field at Mq = 0.5 are shown 
in Figure 5. Four flows are present in this example: 1) the free-stream 

or external flow at Ho = 0.5, 2) the inlet flow for a mass-flow ratio of 
1.0, 3) the fan nozzle flow, and 4) the core nozzle flow. The fan nozzle 
and core nozzle flow have temperature and pressure profiles typical of this 
type of jet engine. 

The details of the fan nozzle flow at a supercritical nozzle pressure 
ratio, are shown in Figure 6, The external flow is at = 0.02 to represent 
a near- static nozzle expansion. The radial shifting of the streamlines due 
to flow curvature is very evident, 

4.2 TWO-DMENSIOMAL IISLET 

The two-dimensional inlet with the ramp adjacent to the aircraft fuse- 
lage (typical of two-dimensional fighter-type inlet with boundary layer 
bleed) is shown in Figure 7, The flow field for this inlet was calculated 
with three flows. The initial free-stream Mach number was M© = 0,8. The 
inlet flow was choked at the inlet throat. 

The pressure distribution on the external surface of the cowl is 
shown in Figure 8. The maximum surface Mach number was 0.998. 

The pressure distributions on the ramp and the upper wall of the 
inlet are plotted in Figure 9. The local Mach numbers exceed unity as 
the ramp turns toward the axial direction in the throat. The walls down- 
stream of the throat were defined as straight in this example, and the 
pressure levels show uniform flow at Mach s 1.0. 

4.3 DATA COMPARISONS - AXISTMMETRIC INLETS 

The NASA inlet No. 8 (NASA 1-85-100), with an internal contraction 
ratio of A m. /A»r^yQ^,^ = 1,093, was selected as the configuration for data 
comparisons. This inlet represents a typical flight-type installation for 
high transonic flight Mach numbers. The geometry consists of a NACA 1 
series external contour with an x/D^ax. = ® OiiLAHnax. = 0.8535 as 

shown in Figure 10. The projected area distributed is shown in Figure 11 • 
This inlet was tested in the 16- foot transonic wind tunnel at NASA- Langley. 
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Figure 5. Short Duct Fan Installation Flow Field 
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Figure 6. Fan Nozzle Flow, Supercritical Nozzle Pressure Ratio 
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Figure 7 . Two-Dimensional (Planar) Inlet 









Figure 11. NASA 1-85-100 No. 8 Inlet Projected Area Distribution 





The Inlet No. 8 was analyzed with the streamtube curvature analysis for 
three Mach numbers and two mass flow ratios. The selected Mach numbers 
were nominally 0,8, 0.85, and 0.90, and the mass flow ratios were nominally 
0.8 and 0.88, The actual values were set to match the measured test results 
for Mach number and mass flow ratio. No mass flows below 0.80 were considered 
because inspection of test data Indicated that flow separation was present 
on the external surface lip. 

The STC analysis Is an inviscld flow prediction. Viscous effects such 
as boundary layer displacement thickness and separation have to be 
accounted for In a separate analysis. A turbulent boundary layer technique 
(SAB) has been incorporated in the transonic analysis (see Addendum). Some 
of the predicted results shown here will Include viscous effects. These 
results will be discussed more fully In the Addendum. 

A typical predicted flow field, after 12 grid refinements, is shown in 
Figure 12 for the Inlet at Mq = 0.92 and a mass flow ratio of 0,88. This is 
not one of the Mach numbers Included above, but is used for demonstrating the 
calculation results. The extra grid refinement in the lip stagnation region 
is evident. 

The comparison of the predicted pressure on the cowl surface with the 
measured pressures is shown in Figure 13. The measured pressures consist 
of three lines of static taps at three circumferential positions (0* , 90® , 
and 180® forward looking aft). Note that the flow is symmetric over the 
nacelle cowl. The comparison plot is arranged so that the highlight diameter 
occurs at an axial distance of zero, and the internal surface of the cowl 
lip is shown as a negative distance. Thus, the surface pressure distribution 
can be represented as a continuous curve. 

The predicted results from STC show that the stagnation point is located 
exactly, and that the pressure distribution is predicted quite accurately. 

The sonic pressure coefficient, Cp*, is indicated, and the comparison shows 
that a shock is present on the inner surface. The mass flow ratio is rela- 
tively high so that there is little acceleration over the external surface. 

The comparisons of predicted pressure distributions from the STC 
analysis with the NASA-Langley test results from the 16- foot wind tunnel 
are shown in Figures 14 through 20. These cover the range of Mach numbers 
and mass flows listed above. 

At Mq = 0.8, the comparisons are generally good for a solution repre- 
senting 700 grid points (solid lines in Figures 14 and 15.) When 1100 grid 
points were used for a mass flow ratio of 0.81, the local oscillations in 
the Inviscld flow were evident. The test data indicate that viscous effects 
on the wall eliminate this pressure fluctuation. Later analysis with the 
viscous analysis predicted local separation (see Addendum). 

At Mq = 0.85 and a mass flow ratio of 0.8819, the correlations between 
measured and predicted still show good agreement (Figure 16), For a mass 
flow ratio of 0,8064 at Mq = 0.85, the effects of compression waves and 
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Figure 12. NASA 1-85-100 No. 8 Inlet Predicted Flow Field 
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Figure 14. NASA 1-85-100 No. 8 Inlet Nacelle Pressure Distribution® Mq = 0.8021 
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Figure 16. NASA 1-85-100 No. 8 Inlet Nacellle Pressure Distribution @ - 0.8510 



viscous Interactions are evident (Figure 17). The Inviscld pressure dis- 
tribution predicted by the STC analysis Indicates compression waves, both 
with 700 grid points and with 900 grid points Including viscous effects. 
Also, local separation is predicted after the first compression wave. From 
the experimental pressure measurements. It Is not evident that the flow 
separation exists over a large region. Once the flow reattaches, there is 
excellent agreement between the measured and predicted pressures. In the 
local region on the cowl lip, there is evidence of viscous shock interac- 
tions, and the Inviscld flow calculated by the STC analysis will not cor- 
relate with the test data. 

The comparisons at Mq = 0.90 indicate that viscous effects are more 
Important. The predicted flow field at a mass flow ratio of 0.885 agrees 
much better when boundary layer displacement thickness Is Included (Figure 
18). At a mass flow ratio of 0.81 (Figure 19), the predicted pressure dis- 
tribution shows several strong compression waves, and separation is indi- 
cated. From the experimental pressure measurements, there Is evidence that 
the flow separation exists over the initial portion of the cowl lip. The 
measurements appear to show a gradual recompression followed by a weaker 
wave at an axial distance of 1.5. The viscous effects are very evident, 
and the inviscld analysis needs to be augmented with a detailed boundary 
layer analysis. 

For a fully subsonic flow, the inviscld analysis by STC predicts the 
wall surface pressures excellently. Figure 20 shows a comparison of Mo = 
0.70 and a mass flow ratio of 0.87. 


The integrated pressure forces, both measured and predicted, are sum- 
marized in Table I. The Integrated pressure drag on the external surface, 
normalized by free- stream dynamic pressure and maximum nacelle area, is 
noted as CDp. For STC, the integration starts at the calculated stagnation 
point and extends over the external surface to the maximum diameter. The 
integrated pressure drag from the NASA-Langley test data is the sum of the 
pressure integrals from the three rows of static taps (0® , 90“ , and 180“ ) 
applied to the complete nacelle (where the 90® row is assumed for the 270® 
row). The equation is as follows: 



P A„ , 
oj ProJ 

%ax. 


( 10 ) 


The total pressure drag on the Inlet is predicted by summing the 
pressure Integral on the stagnation streamline (additive drag) and the 
pressure Integral on the external surface. This total force is again 
normalized by the free- stream dynamic pressure and the maximum area. 
There are no comparable measured drag data, since the force balance in- 
cludes the friction force on the external surface. Also shown in Table 
1 are the number of grid points used in the solution and whether the 
viscous effects (SAB) were included. 
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Figure 17. NASA 1-85-100 No. 8 Inlet Nacelle Pressure Distribution® Mq = 0.8497. 
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Figure 18. NASA 1-85-100 No. 8 Inlet Nacelle Pressure Distribution @ Mq = 0.9007. 
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Figure 19. NASA 1-85-100 No. 8 Inlet Nacelle Pressure Distribution @ Mq = 0.900. 
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Figure 20. NASA 1-85-100 No. 8 Inlet Nacelle Pressure Distribution @ M = 0.6974 










Table 1. Integrated Pressure Forces 


M 

o 


CDp 

Measured 

CDp 

Predicted 

Total 

CDp 

Grid 

Points 

SAB 

0.6974 

0.8715 

-0.023 

-0.034 

-0.008 

715 

No 

0.8021 

0.8754 

-0.025 

-0.023 

0.006 

626 

No 

0.8008 

0.8093 

-0.039 

-0.044 

0.009 

656 

No 



-0.039 

-0.040 

0.010 

1091 

No 

0,8008 

0.8093 

-0.039 

-0.047 

0.004 

997 

Yes 

0.8510 

0.8819 

-0.025 

-0.030 

-0.001 

758 

No 

0.8497 

0.8064 

-0.044 

-0.053 

0.003 

906 

Yes 



-0.044 

-0.059 

-0.028 

717 

No 

0.9007 

0.8852 

-0.0257 

-0.022 

0.009 

700 

Yes 



-0.0257 

-0.031 

-0.001 

928 

Yes 

0.9001 

0.8073 

-0.050 

-0.054 

0.009 

782 

Yes 
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Overall, the inviscid analysis by STC is in general agreement with the 
test measurements, except where viscous interactions are significant* The 
need to include viscous effects is obvious. Further work is necessary to 
properly define the separation bubble and the point of reattachment. 

4.4 AFTERBODY WITH SHOCK 

A nacelle afterbody with a 24® boattail angle was analyzed with the 
Streamtube Curvature Analysis. The gecmietry, shown in Figure 21, represents 
the high subsonic cruise configuration of a nacelle afterbody designed for 
supersonic operation. The afterbody model was tested with a sting^mounted 
forebody at a Mach number of 0.90. A comparison of the experimental results 
and the predicted pressure distribution is shown in Figure 22, 

The predicted pressures indicate a pressure drop or local acceleration 
around the radius onto the 24® boattail, then a sharp compression or shock. 
The location of this compression is a result of the numerical star- switching 
procedures built into the STC analysis, when the velocity changes from 
supersonic to subsonic. The strength of the shock is related to the local 
change in curvature across the orthogonal line defining the star switch. 

The exact Rankine-Hugoniot equations are not included in order to not over- 
constrain the problem. The entropy rise across the shock or compression 
wave on any streamline can be defined, based on the static pressure rise 
defined by the flow field. In this particular case, the upstream Mach 
number was 1,31 and the downstream Mach number was 0.78. This corresponds 
to a normal shock at the body surface. Thus, without including the exact 
shock relations, the predicted compression represents the location and 
strength of a normal shock. 

The experimental results show that the predicted compression wave or 
shock is correctly positioned, but that the boundary layer separated at the 
shock. Once the flow separates, the static pressure is nearly ambient over 
the remainder of the boattail and no effect of the remaining boattail 
geometry or Jet plume compression is evident. Thus, the STC analysis can 
locate the shock position and specify its approximate strength. 
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Figure 22, Afterbody Pressure Distribution 





5.0 DETAILS OF THE NUMERICAL PROCEDURE 


In this section the details of the numerical procedures employed In the 
Streamtube Curvature program are presented. 

5.1 GRID COORDINATE SYSTEM 

When the original grid Is established (as shown In Figure 4) §i and 
§2 coordinates are assigned to each orthogonal line and streamline. These 
values remain attached to the same line throughout the calculation procedure. 
§1 generally has a value of zero at the upstream boundary and Increases by 
8.0 or 16.0 across each region. Similarly, the §2 coordinate Is zero on the 
lower boundary and Is Incremented by 8.0 across each channel. Double stream- 
lines are used to separate the channels and each has the same value of § 2 * 

As the grid Is subdivided during the refinement process, the new lines 
are given coordinate values half way between those on either side. As a 
result, each grid point has a § 1 , §2 coordinate. However, these coordinates 
are for notatlonal and bookkeeping purposes only. §1 and §2 values do not 
enter into the solution of the equations. They are employed in the STC 
Program because, with the conventional counting system (such as the stream- 
line number or the orthogonal line Index) the value associated with a given 
line would be changed when new lines are inserted Into the field, whereas 
the § 1 , §2 values are not. 

The values always increase in the downstream direction, and the §2 
values always increase when one proceeds across the field (to the left after 
facing downstream). Because of the possibility of multiple channels, the 
streamline and orthogonal line Index values are not so ordered. 

Reference will be made to the §1 coordinate in Section 5.3, where the 
§1 value is used to establish the relative spacing between orthogonal lines. 

5.2 CURVATURE OF THE STREAMLINES 
5.2.1 The Beam Fit 

The third step in the calculation procedure as outlined In Section 3 Is 
the determination of the streamline curvatures, angles, and cumulative curvi- 
linear lengths at each grid point. An accurate and rapid method for accom- 
plishing this Is to fit a draftsman's spline or, equivalently, to use the 
formulas which apply to a beam loaded at discrete points. The classical 
relation which Is applicable here Is that the moment, M, varies linearly, or: 


^=^T = bX 

dx2 


( 11 ) 
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In Equation 11, y is the vertical displacement, E is the modulus of 
elasticity, 1 is the cross-sectional moment of inertia, and b is a constant. 
A limitation on Equation 11 is that: 



To ensure that this condition is satisfied, the coordinate system in 
which the equation for the beam centerline is written is rotated as shown 
in Figure 23, As indicated, the curve-fit equations utilize a different 
coordinate frame for each interval. For the interval between point i and 
i+1, the origin is placed at point i and the x-axls passes through point 
i+1. 


The displacement equation for a beam with point loading is a cubic 
polynominal and may be expressed as: 

2 3 t h 

y = b^x + c^x + d^x (for the i internal) 

A more convenient form of the cubic which passes through both points i and 
i^ 1 is: 

i; " (-8^ + - vi (*" - (12) 

where : 

y^ s Slope at X = 0 

y^^ = Slope at X = 6x 

f = x/Ax 

g = (Ax - x)/Ax = 1-f 


The cubic fit to the adjacent Intervals must be matched so that, at the 
Junction points, the angle and curvature are the same. As is illustrated 
in Figure 23, this requires: 


tan"^ [yb(i-l)l 

If 

^b(i-l) 


1 + y 


/2 
b(i-D] 


3 ^ 

|2 1 + y^ 


"a(i 


(13) 

(14) 


Good approximations to Equations 13 and 14, providing Aa is 10® or 
less, are: 


I 
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Figure 23. Curve-Fit Algorithm. 


yb(i) = ya(i+i) + ^a(i+i) 


(15) 


and: 


, 3 /2 

1 + 2 yb(i-l) 


If 


1 


+ 


^ad) 

3 /2 

2 ^a(i) 


(16) 


Substitution of Equation 15 into Equation 12, and the result dif- 
ferentiated gives: 


^ =1^ jj'ad) * \<(UU * ^(id)] <«-«j 


(17) 


At point i, but for the i-1 interval, it follows from Equation 17 


that: 


'^b(i-l) Ax(i-i) 


^>'.(1-1) * ^ |<(i) " 


( 18 ) 


For the same point, but for the interval, we have: 


^a(i) "" 


-4ya<i) - 2 l^add) * ^“d+l)| j 


(19) 


Equations 18 and 19, substituted into Equation 16 and rearranged, give 
the following recurrence equation for ya(i) = • 

^li *^i ^(^li ^3i) *^i + ^3i **(i+l) ^ "^li " '^3i ^°^(i+l) 

(i = 2,3,4 N-1) 

where: 

*11 = h * I ^d)l 


*31 **(1-1) 1^*2 ’’bd-l)) 


(21b) 


Note that b^ is used in place of ya(i) simplify the notation. 

Notice that because of the ya and dependency, the values of A cannot 
be evaluated directly; Instead, they are determined by iteration. However^ 
because y^ is always very small, two passes (one correction pass) are 
generally sufficient for their determination. 
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Given and together with the boundary condition equations pre- 

sented below, the set of equations given by Equation 20 yields a tridiagonal 
matrix equation which is solved by standard procedures, y^ is then evaluated 
using Equation 15. Likewise, the ordinate at any intermediate point can be 
computed by Equation 12; and Equation 16, together with Equation 19, gives 
the curvature at each point. 

The curvilinear distance over the interval is given by the expression: 

i+1 1 

f ds = Axj^ f 

The first two terms of the binominal expansion are: 

^1 + y'j^ = 1 + iy' + . . . • (23) 

Using Equations 12 and 23 in Equation 22 yields: 

®(i+l) " ^ ^^ayb + yh^)A^\ <24) 

This equation is employed to calculate the cumulative curvilinear dis- 
tance along the streamlines. 

5.2.2 Beam End Condition Options 

Three different end options are available with the beam curve fit: 

• The angle may be specified. 

• The curvature may be specified. 

• The ratio of y”* (the rate of change of curvature) at the end point 
to the value at the next- to- end point may be specified. 

As a ground rule, the second option is used for fitting the streamlines. 
For those streamlines which extend to the flow inlet or flow exit boundary, 
the end curvature is taken as zero. (A constant value of curvature, dif- 
ferent from zero, can also be enforced by user input.) If a streamline 
terminates within the field, the end curvature is interpolated from the 
curvature of the streamlines above and below. 

The equations for the three end options, derived from the formulas 
listed in Section 5.2.1, for the first end are given below: 

a) Specified Angle: 
bj = tan ^0^ " 




df 


( 22 ) 
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b) Specified Curvature: 

4bj^ + 2b2 = "2^0.2 + ^ 2 ^^d)] 

c) F specified where y^'” = F 72*” 

+[(AX2)2 - F(AXi)2]b2 - F(AXi)2b3 

= F Att3(Axp2 _ Aa2<^2)^ 

where : 

01 = Specified angle 

2 2 

Cl = Specified curvature (= - d y/dx ) 

ai = Angle of the chord between points 1 and 2. 

The parameter F can be given the following interpretation. When F = 0, 
the curvature in the end interval is constant. Hence, the end interval 
polynominal in this case reduces to a parabola. When F = 1, the third 
derivative in the first interval is equal to the third derivative in the 
second interval. In this case, since angle, curvature, and y*” are all 
continuous at the second point, the same cubic equation spans the first 
two intervals (at least to a good approximation for small Aa 2 ). 

Similar equations to Equations 25, 26, and 27 can be written for the 
downstream (or second) end of the beam fit. However, they are omitted here 
for brevity, 

5,2,3 The Stagnation Point End Condition 

The curvatures for points on a body surface are computed from the input 
geometry data, as will be discussed in another section. The curvature at a 
trailing edge point is taken to be the same as the body surface curvature; 
although, in reality, there may be a weak singularity at this point because 
of a finite \yedge angle. Just downstream of the trailing edge, if the flow 
is subsonic, the curvatures are computed by the beam formulas where, at the 
trailing edge point, the third option (F = 0) is used for the upstream end 
condition. 

The leading edge is handled in a slightly more complex way. It is a 
requirement in the STC Program that the leading edge be rounded and that a 
complete numerical description of the leading edge shape be supplied. At 
the stagnation point, then, the streamlines are required to turn a 90® 
corner. Two coincident stagnation streamlines are employed. One turns up 
and goes over the body; the other turns down and goes below. The location 
of the stagnation point is found iteratively. As shown on Figure 24, the 


( 26 ) 


(27) 
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Figure 24 . Stagnation Point Iteration. 



point is moved along the contour so that the Intersection angle with the body 
surface is 90" . The streamline angle at the stagnation point is found by 
utilizing the spline fit with the third option (F = 1) for the downstream 
end condition. 

5.2.4 Backward Curve Fits for Supersonic Regions 

When the flow at a particular grid point is supersonic, the "complete" 
beam fit (or central difference) formula is not appropriate for the evalua- 
tion of curvature. To be consistent with the physical character of super- 
sonic flow, the curvature at point 1 cannot be influenced by any points down- 
stream of that point, such as the i+1 or i+2 points. 

Therefore, to evaluate curvature at supersonic points, the "beam" is 
fitted to only 3 or 4 points to obtain a curvature at the last of these 
points. The beam-fit is moved along (by dropping a point at the upstream 
end and picking one up downstream) to calculate the curvature at sequential 
points on the streamline. Used in this way the beeun fomulas essentially 
represent either a parabolic (3-point) or a cubic (4-point) curve fit. The 
cubic fit is obtained when F at both ends, for the 4- point fit, is set to 
unity. Actually, more favorable agreement with theory (in the case of a 
two-dimensional Prandtl-Meyer turn) is obtained when F is set to about 0.75, 
and these are the standard values of F present in the STC program for the 
4-point formula. F is automatically set to zero when the 3-point option is 
selected. 

For pure supersonic flows, the 4- point curve fit (which is second 
order accurate) gives much better results than does the 3- point (first 
order accurate) formula. But for mixed flows, it is found that the 4-point 
formula generally leads to divergence and, therefore, the 3- point formula 
is always used for transonic cases. 

5.2.5 The Evaluation of Curvature Very Close to a Sonic Line 

In the STC method, the curvature must be known at the sonic (or near 
sonic) points Just as at any other i>olnt in the field for use in the cross- 
stream momentum equation. Equation 7. This is in contrast to the velocity 
potential method of Murman and Cole where the second streamwise derivative 
of <t> is insignificant because the coefficient (1-M^) is zero. Therefore, 
the following procedure is used to evaluate the streamline curvature at 
points close to a sonic line. 

The point at which the sonic line crosses each orthogonal is identified 
by referring to the value of velocity (or, in reality, the coefficient B) 
from the previous iteration. Then, for some distance above and below the 
sonic line, the curvatures are taken to vary linearly between the values cal- 
culated according to Sections 5.2.1 to 5.2.4. This is illustrated in 
Figure 25. The extent of this linearly assumed variation is controlled 
by an input variable (TSIC). 
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5 . 2.6 Special Interior Points at Orthogonal Ends 


As mentioned earlier, it is not required that a newly inserted orthos 
onal line span the field if the numerical resolution near the boundaries 
is already satisfactory. In this case the orthogonal line will terminate 
within the field on a streamline. These are considered special points and 
they are not used in the curve fits described in Sections 5.2,1 and 5.2.4. 
Instead, after the curve fit is obtained, the position, angle, and curvature 
at the special points are interpolated. 

Another place where this procedure is used is in the positioning of 
points adjacent to a stagnation point, as indicated by the solid circle in 
Figure 24. In this case the interpolation procedure replaces the integra- 
tion of the momentum and continuity equations in the interval containing 
the singularity at the stagnation point. 

5.3 POSITIONING THE ORTHOGONALS 


The cross- stream momentum and continuity equations are written in a 
direction normal to the streamlines. Hence, before these equations are ap- 
plied, it is necessary to move the grid points along the streamlines to ob- 
tain orthogonality. This is easily accomplished as follows (refer to Figure 
26) . In each region a boundary or dividing streamline is chosen as a ’’control** 
streamline. Along this line the spacing between orthogonals is chosen pro- 
portional to A§i. To correct the nonorthogonality, the points on each orthog- 
onal are first fitted with a spline (using the equations of Section 5,2 and 
the end condition F = 0) to obtain the angles 02. The angle deviation from 
the streamline normals is then integrated with respect to the cross- stream 
distance n, from the control streamline to the point in question, to obtain 
the relative movement Ds, The points are then moved to the new positions by 
utilizing the streamline curve- fit equations. The coordinates, 0i angles, 
curvatures, and cumulative s-distances are modified as appropriate. The 
constant of integration, Ajj3, in the integral for Dg positions the orthog- 
onal on a control streamline so that a reference position is maintained. 

5 . 4 FAR-FIELD SOLUTION 


The boundary condition on the external flow is that the velocity approach 
the undisturbed velocity, V®, and the flow angularity approach zero as the 
spacial coordinates approach infinity. To economically meet this condition, 
the field is divided into an ’’inner** and an ’’outer” region, illustrated in 
Figure 27. The inner region is the region near the body where flow distur- 
bances are large and/or the typical nonlinear transonic effects are encoun- 
tered. Flow properties in this inner region are calculated by the streamtube 
curvature technique which uses the full nonlinear equations of motion. 

The outer region is the region extending from the outside edge of the 
STC integration domain to infinity. In this region, an asymptotic form of 
the equations of motion is applied. These asymptotic equations are solved 
analytically. On the interface boundary between the two regions, it is re- 
quired that the velocities, as calculated in the two separate regions, must 
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Figure 26 . General Positioning of the Orthogonals and Point Movement for Obtaining 
Lines Normal to the Streamlines. 




match in both magnitude and direction. Consequently both the velocity on the 
boundary and the shape of the boundary are matched. This manner of applying 
the boundary condition is analogous to the familiar inner-outer expansion 
method of asymptotic theory (ref. 6). 


To achieve a matched solution on the interface boundary, there are 
really two questions to be answered. First, for a particular interface 
streamline shape, what is the velocity distribution? The determination of 
this velocity distribution is the fifth step in the calculation procedure 
outlined in Section 3,0, and these velocities are used as boundary conditions 
on the momentum equation in the sixth step. The second question is: if the 

interface streamline position is altered slightly (by moving the points out- 
ward or inward), what will be the change in velocity at any given point? 

This is required in the matrix formulation of the streamline correction equa- 
tion, Section 5,6,4, In this section only the first question is considered. 


In the present procedure, the far- field region is approximated by a 
linear differential equation, the solution to which can be computed very 
rapidly. The linear formulation, however, also requires that the distur- 
bance of 3^ be small in the outer region compared to the free- stream value, 
3», where 3^ = (l-M^), Thus, the far-field boundary must be sufficiently 
far from the body to satisfy this condition, (Alternate far-field solution 
procedures which allow a relatively large 3£ disturbance could, perhaps, be 
employed as mentioned at the end of this section,) 


The calculation method is a simplified version of the Douglas-Neumann 
procedure which was developed by Smith and co-workers (refs. 7 and 8) at the 
Douglas Aircraft Company. Two separate versions of the analysis have been 
formulated, one for axisymmetric flow, the other for two-dimensional flow. 
Both of these analyses are similar in nature, but the mathematics of the two- 
dimensional solution is much simpler. In fact, the two-dimensional analysis 
degenerates to a simple thin airfoil calculation. Because it is simpler, the 
two-dimensional analysis will be presented first. 

Since the **body’* of interest in the outer region really represents the 
outer streamline of the STC numerical integration domain, it is reasonable 
to expect that its ’^thickness** (i.e., the vertical displacement of the inters 
face streamline frcati its undisturbed state) will be small. Further, all 
angles on the *’body” surface (i.e., flow angles on the interface streamline) 
are also expected to be small. The assumption of small deflections and 
angles and the assumption of small transonic effects allows us to use the 
classical small perturbation equations of subsonic flow theory to describe 
the outer flow field. Thus, the problem of finding the velocity on a surface 
which is only slightly perturbed from a straight line is just the thin air- 
foil problem. 

As discussed in any aerodynamics text (for example, see ref, 9, Chapter 
7), the velocity potential for the flow past a thin airfoil is governed by 
the equation: 
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0 


(28) 


p2 


which can be converted to Laplace's equation, y2 | _ q, by a Prandt 1-Glauert 
transformation. The solution, for the thin airfoil approximation. Is obtained 
by putting a series of mass sources on the original unperturbed straight line 
and matching the velocity tangency condition there rather than on the actual 
perturbed boundary which describes the body shape. We can then write the 
solution for the velocity potential, §, and the perturbation velocities u-ito 
and V as: 


$(x,y) 



(29) 


U-Uco = 


1! 

5x 


=J- 

PtT j._ri 


g(g) (x-C) 


LE (x-|)^ + 8^(y-'n)^ 


dC 


- - i ( 

" Sy “ n J, 


TE 


LE 


ct(C) 8(y-n) 


(30) 


(31) 


where §, r| Is the location of the source point and where x, y is the field 
point at which the potential (or velocity) Is being calculated. 


In the case of the far field, all mass sources are placed on the straight 
line which corresponds to the undisturbed Interface boundary. We are in- 
terested in finding the velocity on this same line. Thus all calculations are 
made on the line T| = ys. (The subscript, s, refers to the location of the 
undisturbed interface streamline.) The velocity distribution on the body 
can be found by evaluating Equation 30 at y = y^, and r| = yg. Thus, we have: 

u(*,„)-u. = ( 32 ) 


We now subdivide the interface streamline into a series of N equally spaced 
intervals. Thus the integral In Equation 32 beccxnes the summation of N 
Integrals: 


1 N xj+^ 

u(x,y,)-u. = -g; ^ c(xj) f 


-J-i 


(33a) 


Performing the indicated integration and writing the field point x as the 
subscripted point Xj^, we have: 
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(33b) 


u(xi,yg)-Uoo 


N 


o(K^) In 


J 



The sum of integrals in Equation 33a includes one integral in which the 
integrand is singular, namely that interval in which lies. By considering 
the Cauchy principal value of the integral in this interval, it can he shown 
that the singular interval's contribution to the velocity is zero. For this 
reason, the summation in Equation 33b indicates that the interval j=i is to 
be discarded when calculating the velocity. 

A shorthand notation for some of the terms in Equation 33b is now intro- 
duced. Let: 



so that Equation 33b becomes: 

- «00 = Xij CTj (35) 

where u^ and qj represent u(xi,ys) and a(xj), respectively. 

Determination of the source density is simple in the two-dimensional 
case. The density is given by: 

aj = 2Vj (36) 

where v is the vertical component of velocity on the body. This relationship 
can be obtained frcxn Equation 31 by a careful consideration of the limiting 
process as y approaches the location of the source, r|. Again, the derivation 
can be found in standard aerodynamic texts. By using the small perturbation 
assumptions: 


V 



(37) 


Finally, for reasons which will become apparent in the axisymmetric analysis. 
Equations 36 and 37 are combined to define the matrix such that: 





where is a diagonal matrix. 


( 38 ) 
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We now combine Equations 35 and 38 to obtain a direct equation for the 
u-velocity component: 



(39a) 


or: 


sj 


where , 


Zij = Xik 



(39b) 


(40) 


Of course in the thin airfoil approximation where (like Y^j) is a 

diagonal matrix, Z^j is identical to except for a constant. 

It is interesting to note that if the sources were placed on the body 
surface instead of on the line y = y^, the matrix would be a completely 
dense matrix so that finding its inverse would not be trivial. However, for 
very thin bodies, the off-diagonal terms are very small; i.e., the matrix 
is strongly diagonally dominant; and, in the limit of a very thin body, 
will approach a diagonal matrix. (In the axisymmetric analysis, the 
matrix Yij will also be a completely dense matrix.) Placing the source dis- 
tribution on the body surface, makes the two-dimensional analysis identical 
with the Douglas-Neumann analysis (ref. 7). 

During the iteraction procedure for the Streamtube Curvature solution, 
the matrix Z^j is calculated only once. Thus, even though the ’’body** (i.e., 
the interface streamline) changes shape, the matrix Z^j remains unaffected. 
This is because the thin airfoil approximation allows sources and the tan- 
gency condition to be placed on the straight line y = yg, not on the per- 
turbed interface streamline. The entire effect of the body geometry enters 
through the vector, (dy/dx)gj, which represents the slopes of the body 
surface. 

Before the calculation procedure can be implemented, the upstream and 
downstream limits of the integration must be specified. We have chosen to 
place the ’’leading** and ’’trailing” edges of the interface streamline upstream 
and downstream of the numerical (inner) integration boundaries, respectively. 
This spanwise lengthening of the outer domain over the inner domain is done 
to ensure that the velocities which are calculated by the outer analytical 
solution are reasonable and well-behaved at the streamwise extremities of the 
STC (numerical) integration domain. To allow the streamwise extension of the 
interface streamline, a quadratic addition is fitted to both ends of the STC- 
calculated interface streamline. These quadratics are defined by requiring 
that they have both zero deflection and zero slope at their outer ends, and 
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that they join continously (but allowing discontinuities in slope) with the 
ends of the STC streamline. 

The analysis of the axisymmetric version of the analytic outer region 
solution is identical in philosophy to that of the two-dimensional solution. 
The generation of the matrices ^ij» ^i j * however, is more compli- 

cated. The details of the axisymmetric analysis have been given by Smith and 
Pierce (ref. 7). Consequently, only the results will be included here. 

In axisymmetric flow, the interface streamline cannot be thought of as a 
thin airfoil, and it certainly cannot be considered to be a slender body. 
Instead, the interface streamline appears as an outward perturbation on what 
would otherwise be a right circular cylinder. As in two-dimensional flow, we 
assume that a series of mass sources can be placed on the straight surface of 
the cylinder itself rather than exactly on the interface streamlines. These 
sources must not be simple mass sources, but rather must be ’Ving” sources. 
The potential field which is induced by a ring source of radius a, at the 
axial location Q, is: 


d$(Z,r) = 2ac7d<; 


^ d0 

^o [(z-Q^ + r^+a^ - 2ar cos 9 ]^ 


(41) 


By placing a continuous distribution of these sources on the cylinder r = a, 
and then breaking it up into a number of small segments as in the two- 
dimensional case, we obtain the following equations for the matrices X^j, 

Y. . (see refs. 7 and 8): 


j+i 


Xij = -4a 


Xij = 0 




E(m) dC 


j-i [4a2 + (Zi_^)2]^ 


and: 


J+i 


r = -2 C K(m) - E(m) 

^ [4a2 + (zi“p2]i 


dC 


D / D\ 

Yij =-2 tt+ 2D in--- ( 1 + in-) 


(i5^j) 

(l=j) 

(i=j) 


(42) 


(43) 


where : 


m = 


4a2 

4a'2 + (Zj^-C)^ ’ 


z . 1 - z . 1 

1+2 i-i 
2a 


D = 


and where K(m) and E(m) are the complete elliptic integrals of the first and 
second kind, respectively. 
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As in the two-dimensional theory, we now calculate the u component of 
the velocity in the outer region from Equation 39b; 



(repeat of Equation 39b) 


where, as before: 


*ij = *ik KiS 


(44) 


Note in this axisymmetric case, that is not a diagonal matrix, and that 
y^j must be found from a numerical matrix Inversion algorithm. However, the 
only geometrical terms which enter the Y^j matrix are the undisturbed radius 
of the outer streamline and the total axial length of the outer region. 

Again, as in the two-dimensional theory, the effects of the displacement of 
the interface streamline enter only through the vectors (dy/dx)gj. This is 
a direct result of placing the ring sources at the constant radius a, rather 
than exactly on the interface streamline itself. Consequently, the 
matrix can again be calculated once for an entire STC iteration history, 
despite the fact that the shape of the interface streamline is changing. 

It is helpful to compare the approximations in the above formulation to 
those of Murman and Cole (refs. 1 and 10), Their formulation is similar to 
the one described above, except that their solution domain includes both the 
near and far fields. They include the nonlinear term by means of a series of 
transonic sources distributed throughout the numerical integration domain. 

But in the region which corresponds to the far field, the transonic source 
terms are neglected. (In fact, they cannot be included because the infor- 
mation required to compute the distributed source strength is not available.) 
Therefore, their method, like the method presented here, neglects the non- 
linear effects in the outer region. 

Despite these justifications, some inclusion of the nonlinear effects 
in the far-field solution would be advantageous. The local linearization 
method of Sprieter (refs, 11 and 12) would allow these effects to be in- 
cluded in an approximate manner, and could be included in future program 
development activity. 

Typical results of the far-field matched solution are now illustrated. 

As anticipated, the proximity effect of the numeric field outer boundary 
is small when the analytic far-field solution is employed. Figure 28 illus- 
trates the variation of peak Mach number on a body of revolution as the 
position of the outer boundary is changed. Notice that, for this test case, 
the induced Mach number error is less than 0,01, even when the outer boundary 
is reduced to twice the body radius (Roster “ 1*0) • 

Also illustrated in Figure 28 are the solutions obtained by utilizing 
hard wall and constant pressure boundary conditions as a function of Router 
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Peak Mach Number Vs. 1/R 
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Figure 28. Peak Mach Number for Far-Field, Constant Pressure, and Solid Wall Boundaries 







By comparison, the utilization of the far-field solution allows a significant 
reduction of boundary radius, for the same accuracy, yielding an efficient 
computational procedure. 

5.5 INTEGRATION OF THE CONTINUITY AND MCMilENTUM EQUATIONS 

In this section the numerical procedure for integrating the full non- 
isentropic and variable energy form of the equations of motion is developed. 
Only a portion of this generality is required, since external flows are 
isoenergetic and only slight entropy variations arise downstream of shock 
waves. For completeness, however, the formulation for general flow proper^ 
ties is retained. 

The ’’heart” of the STC method is the integration of the momentum and 
continuity equations, which are repeated here: 

Momentum: 


1 1(V^ ^ «CV2 + is - 

2 bn “ an 


(2b) 


Continuity: 



( 1 ) 


Equation 2b is the basic Crocco form of the momentum equation containing the 
variables H and S (enthalpy and entropy). These variables may be replaced 
by total temperature and total pressure for greater engineering convenience. 
If it is also assumed that the specific heat Is constant, then the following 
equation is equivalent to Equation 2b: 


1 d(V^) ^ . i Vf ^ ^PT 

2 dn ~ 2 dn 8n 


(45) 


The above momentum equation is integrated by parts. First, we rewrite 
Equation 45 as: 



dP 

(-2Cdn + dlnT^) + R^T 


(46) 


The formula for integrating by parts is: 

dv = - — du + — d(uv) (47) 

u u 

By comparing Equation 47 with Equation 46, it follows that: 

V = ^ (48a) 
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(48b) 


~ = Xdn - dlnT™, 
u ^ 


and, 


- d(uv) = R„T 
u & 



(48c) 


Since the total temperature and streamline curvature are 
Equation 48b can be integrated from the first to the k'th 
obtain: 

^ “■> 

“k = e (ui = 1) 


assumed known, 
streamline to 


(49) 


Equation 48c is then integrated to yield an expression for the product 
of u and velocity squared. 


1 

2 


u 

k k 


1 

2 


u 


k+1 k+1 



u HgT 



(50) 


Equation 50 represents the integration across one streamtube bounded 
by streamlines k and k+1. The integration is performed from the outside 
toward the center because, frequently, the velocity is known on the outside 
boundary. The finite-difference form of Equation 50, employed in the com- 
puter program, is: 


‘k^k 


k+1 k+1 


U’-V' u, ,V. , ■s/'V*k^ 



(51) 


Note that this expression is explicit when no total pressure gradient exists; 
otherwise the expression (Equation 51) is implicit because the static tempera- 
ture which appears on the right-hand side is a direct function of velocity. 



However, the implicit nature is very weak up to transonic speeds, and a 
simple successive approximation procedure is utilized to update the right-hand 
side until the computed fractional velocity change (on each streamline) is 
less than 1 X 10"5, 
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The numerical integration of the momentum equation, then, is a two-step 
process: first the integration of Equation 49 is performed followed by the 

integration of Equation 50, In Equation 49, the integral of curvature is 
evaluated by fitting a second-order polynominal in each interval, (The 
second derivative of that polyncmiinal is established by a least square fit 
to the two nearby points,) Implied in the integration is the fact that total 
temperature and total pressure are known as a function of the streamline index 
k and, hence, as a function of the cross-stream distance n^. This is in slight 
contradiction to the streamwise momentum and energy equations. Equations 3 and 
4 of Section 2,0, since total properties are correctly related only to the cumu- 
lative flow rate (Y = W), However, this slight error is automatically corrected 
in the later stages of the iteration when the assviraed streamline positions are 
effectively coincident with the true streamlines. 

As noted above, the integration starts from the outer boundary and pro- 
ceeds inward. For external flows, the velocity on the outer boundary is ob- 
tained from the far-field equations of Section 5,4, For internal channels, 
this outer boundary velocity is found iteratively as will be discussed below. 
Velocity changes across slip lines are obtained by Equations 8a, 8b, and 8c 
as discussed previously. 

Following the momentum equation, the continuity equation is integrated 
by the algorithm: 


^k+1 ” ^k ^ 


k+1 k 




(53a) 


The average value of mass flow per unit area, in the denominator of Equation 


may be approximated two ways: 


<pv> = 4 |(pV)^ . (pv)^^ J 

(53b) 

<pv> = [(pV)^ 

(53c) 


If the variation of (pV) between streamlines is small, then the two expres- 
sions give similar results. Although Equation 53b is somewhat faster to 
compute. Equation 53c has been found more reliable for some special cases. 
Consequently, an approximation to Equation 53c is employed in the computer 
code. 


If the cumulative cross-stream areas are being computed for an internal 
station, it is required that the last area, A|{, equal the geometric area of 
the passage at that location. To accomplish this, the value of used as an 
initial condition in the momentum Integration is varied. Figure 29 shows a 
typical variation of Area, Ajj, with velocity, Vjj. Obviously, there are two 
solutions — one for the subsonic branch and another for the supersonic 
branch. Although the choice of branch may be controlled by user input, the 
subsonic branch solution is always employed with an inlet/nacelle configuration. 
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Figure 29. 


Illustration of Method for Choosing Outer Boundary 
Velocity in a Confined Passage. 
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A third possibility is that the geometric passage area is below the 
minimum computed area, Aj^. In this case, the flow is said to be choked, 
and the flow is adjusted so that the minimum calculated area will be exactly 
equal to the geometric area. This logic is utilized, for example, in the 
throat of the nozzle discharge passage. 

5.5.1 Stagnation Points 

If a sharp concave corner is encountered, the flow velocity is known to 
be zero. Such a situation cannot easily be handled with the equations as 
formulated above. Hence, the integration of both the momentum and continuity 
equations is omitted in the interval adjacent to the body, and the first point 
away from the body is interpolated as indicated in Section 5.2.6. 

To replace the omitted equations at the leading edge stagnation point, 
the derived condition that the stagnation streamline intersect the body 
at a 90® angle is utilized as discussed in Section 5.2,3. The streamline 
angle at the stagnation point is, of course, double valued. Consequently, 
the orthogonal line which passes through the stagnation point is made per- 
pendicular to the average of the two angles, a requirement which is derived 
from potential theory. 

5.5.2 Wakes From Blunt Trailing Edges 

Because it is much easier to obtain a valid numerical solution if the 
flow streamlines are smooth and the curvatures are not excessive, a dead 
region is allowed to exist behind a trailing edge which has thickness. As 
shown in Figure 30, the thickness of the *’dead*’ region is gradually reduced 
to zero as one proceeds downstream. The derivative, db/ds, has a nominal 
value of 0.1. 

To include the wake displacement effect, the wake area is added to the 
right-hand side of Equation 53a if a channel boundary streamline is crossed. 

In this way the cumulative streamtube flow area includes the wake displace- 
ments. No correction is required in the momentum equation because pressure 
continuity is always enforced across a slip line. 

5.6 STREAMLINE CORRECTION EQUATION 


5,6.1 General Formulation 

In the STC calculation procedure, we start with an estimate of the 
streamline positions, that is, a set of z^,r^. For the second iteration 
cycle, ZqjTq will be the coordinates determined by the first iteration, and 
so forth. These "assumed" streamline coordinates are used to ccmipute the 
cumulative width of the streamtubes (Uq) , curvature, velocity, and density 
and, then, a second set of cumulative streamtube widths (n^). This is 
illustrated in the following listings: 
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Given Operation 

z^,Tq Streamline Curve Fitting 

z^,r^ Orthogonal Line Curve Fitting 

C^,nQ Momentum Equation 

Vq Energy, Entropy, and Eq. of State 

Po>^o Continuity Equation 


Compute 

Co 

no,A^ 

Vo 

Po 

Ax>“x 


If nx equals Uq, the solution is converged. If n^ is different from no, then 
the streamline positions must be adjusted so that the difference between nx 
and no is reduced. The amount by which the streamlines are to be adjusted 
is denoted by 6n, the equations for which are formulated in the present 
section. 


The calculation formulas to this point have utilized the full nonlinear 
equations of motion. Now we employ an approximate set of linearized equations 
which will provide for the computation of 6n; and, through repeated appli- 
cation, bring the discrepancy (nj^n^) to some small neglectable value. The 
final coverged flow field solution, of course, will be the solution to the 
nonlinear equations (l.e.. Equations 1 through 4). 

The continuity equation for one streamtube may be expressed as: 


^k+1 - = 


k+1 

J pVdA ^ <pV>^^j P 


k+1 



where , 

A = J 2 TT r dn (axlsymmetric) 

A = n (plane) 


<pv> 


k+i 


average flow per unit area for the streamtube 
bounded by the k and k+1 streamlines 


(54) 


For two adjacent streamtubes: 


A 
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k+1 

k+1 
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k 
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k 
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- A 


k-1 


- W 


k-1 



k-i 


(55) 


(56) 
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The previous two equations are subtracted to obtain the second difference for 
A, defined as: 


AW” 





k-i 


w, 


k-1 


(57) 


^ ‘ ■ <p''> k-i 


(38) 


Variables A, p, and V (i.e., those without subscripts) will be used in this 
section to represent the ’’corrected*’ (or true solution) values. Thus, 

Equation 58 represents the ’’correct” solution. 

An equivalent expression can be written to represent the continuity 
equation which was used to predict the x- subscripted areas from the velocities 
and densities associated with the assumed streamlines. That is; 



A W 


1 


1 

^^o^o^k-2 


(59) 


Equations 58 and 59 are now substituted into the following identity: 


A^(A-An) ^ A^x a2a = A2(A^-Ah) 

AW AW " AW AW 


(60) 


to obtain: 


A2(A-Aq) r 1 





<PoVk-i 




AW 


(61a) 


The difference between the "correct" value and the value associated 
with the assumed streamlines is denoted by 6 ( ) . In particular: 

6A = A-Aq 

6n = n-nQ 

6V = V-Vo 

K^pv) " 

6C = C-Cq 


64 



Then Equation 61a can be rewritten; 


AW 



A^(Ax~Aq) 

AW 


(61b) 


This equation represents a difference correction form of the continuity 
equation. 

To eliminate the pV terms in Equation 61b, the momentum equation is 
introduced. In regions where the rotationality is zero, the momentum equa- 
tion is: 

~ -Cdn (62) 

As illustrated in Figure 31, the outer boundary of the field is denoted as 
the N'th streamline; and, it is assumed that, between the streamline and 
the streamline, there is a slip line with index j. Integration of 
Equation 62 from the to the jth streamline and then from the to the 
Nth streamline yields: 

j N N 

InVN - InVja + InVjjj - InVj^ = -J Cdn - J Cdn = -J Cdn (63) 

k j k 


where Vj|^ and are the velocities just below and above the slip line. For 
the assumed curvatures, Cq, the momentum equation is: 

N 

lnV„„ - lnV„ * 1BV„ - lnV„^ = -J C„dn (64) 

k 


Equation 63 is subtracted from Equation 64, and it is assumed that the 
lengths of the orthogonals from k to N are not appreciably different. 




(C-C )dn + ln(^) - 

V^o/jb 



(65) 


Since V = Vq + 


6V we use: 





as a first-order linear approximation, and Equation 65 can be written as 
follows: 



65 



Figure 31. Streamline Subscript Notation 



or. 


{¥) - (€) “ f 

Voile Vo/^ 


( 66 ) 


The parameter aj is introduced to simplify the equations, where: 

- ■ tel - m 


(67) 


'Ja 


Equation 66 is the difference correction form of the momentum equation. 
With it we can approximate the velocity corrections for the (k+ 2 ) and 
(k-|) streamtubes as follows: 



For any given streamline the flow is isentropic, and it can be shown that: 


d 



^ 

pV V 


(69) 


where - )--m2. 


A discrete linear approximation to Equation 69 is: 


-6 



k+^ 



and: 


-6 



k- 



(70a) 


(70b) 
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Equation 70a is subtracted from Equation 70b, and Equation 68 and 68b are 
substituted into the result. This yields the following expression for the 
pV-terms in the continuity difference correction equation: 



Rearranging the above, it follows that: 



(71) 


(72a) 




(72b) 


The term aj, which represents the difference between the velocity 
correction across a streamline, is now to be evaluated. The definition 
of a j , given by Equation 67, is repeated: 



The subscripts "b" and "a" are used to indicate the value of velocity 
just below (b) and above (a) the slip line along streamline j; the subscript 
"o" has been dropped: 



The equations which apply at the slip line are: 


•ivj = c I 

2 b p< 


Tb 


- T 


sb , 


(74a) 




Ta 



(74b) 
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(75a) 


Tsb = Tjb 


T - T 
*sa ” Ta 


&) 

fe) 


Izi 

Y 


111 

Y 


These equations are differentiated to obtain: 

ii-i -1 

Y 


- Vb dVb = -dTgb = 

P V‘Tb> 




dP 


Tb 


Cp dVa = -dTsa = 


fc) 


Y 


dP 


Ta 


Dividing one by the other: 


Vk dV, T^^ 
P b Tb 

V dV “ T^ 
a a Ta 


fe) 


111 

Y 


An approximation to Equation 76 is: 


111 

V 6V^ T jv \ Y 

a b a Tb / Ta I 

V. 6V “ v2 T„ \P„^/ 

b a Ta \ Tb/ 

and Equation 73 can now be written as: 



V2 T„^ i 

111 1 

/6V \ 1 

\ Y 

1 ^1 

a Tb / 

Ta\ 

1 ® ' J 1 

t ’’Ta' 

L b 

.^Tb/ J 


Just above the slip line, the velocity variation is given by the 
of the momentum equation (see Equation 66): 


(75b) 


(76) 

(77) 

(78) 
integral 


69 





6C dn 


So, oj becomes: 


Oj = 


Y 


1 (!s\ ' -1 M . c 

^Ta VTb/ Jj L\ °'N 

b 


6c dn 


= B 


3j 


(€) „ ^ “ -”] 


where: 


B3J = 


Pm 

a Ta 


Y 


P 

b Tb 


Y 


-1 


Ta 


Tb 


Equation 81 is now substituted into Equation 71: 




N 


B 2 B 3 j/^C dn 

And this result is substituted into Equation 61b: 

- Bi 6C^ + SC dn + Ba(l + •» 


6C dn 


A2(Ay-Ap) 

" AW 
where, again: 


Bi 


2 [(poVo)k+^ Nk+1 


" X) (poVo)k_| " "°k-l^ 


(79) 


(80) 


(81) 


(82) 
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Equation 82 is the desired correction equation. The coefficients 
B2) and B^. are all based on flow properties calculated from the assumed 
curvatures. They are known quantities. The unknown quantities are 6A, 

6C, and the velocity variation on the outer boundary 6V|^. Each of these 
quantities must be expressed in terms of 6n, 

Remember that for axisymmetric flow: 

6A = 2tt r^ 6n (83a) 

and for plane flow: 

6A = 6n (83b) 


Thus, for axisymmetric flow, the second difference of 6A is: 


A^(6A) _ / ^k+l 

AW Vl‘\ 


** ^k-1 ^"k-1 \ 
\ ~ \-l ' 


(84a) 


and for planar flow: 


A^(6A) _ ^*Hc+l ~ ^"k _ ~ ^“k-1 


AW 


W, , “ W, 
k+1 k 


k k-1 


(84b) 


The curvature correction is: 
D^(6nk) 


-6C fn 
k ~ 


Ds2 


where s is the curvilinear distance measured along the k^^ streamline. For 
any selected curve fit (or finite difference second-derivative formula), 
the Influence coefficients can be represented as follows: 
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( 85 ) 


D^( 6 nk) 

Ds^ 


% 

G]^6ni^k + Cr26n2^k + G36n3^k + 6^604, k + G 5 &n 5 ^k 

+ G66ri6,k 


Here the numbering scheme shown in Figure 31 has been used. Generally the 
influence coefficients, G, are computed by assuming a spline curve repre- 
sentation. If the velocity is subsonic, a central point difference 
equation is used; and, if the velocity is supersonic, a backward difference 
parabola equation is used. 

With Equations 84 and 85 substituted into Equation 82, there 
results a set of equations for the variables where i is an orthogonal 

line index and k is the streamline index. These* equations, together with 
the boundary condition equations presented below, form a solvable set of 
linear simultaneous equations. 

However, this resulting equation is quite cumbersome. The integral 
terms which appear in Equation 82 lead to a dense coefficient matrix and 
a computer solution is impractical. Therefore, the terms in Equation 82 
which have coefficients B 2 and B 3 are neglected, and the following approxi- 
mate form of the correction equation is utilized: 

= AMS--...AO). ( 86 a) 

AW ^ “1 Ds^ AW 


This form is equivalent to the differential form of the correction equation 
derived in Appendix A, For most flow conditions Equation 86 a produces 
rapid convergence. Hence, the omitted terms seem unnecessary for these 
cases. In a few cases the direct use of Equation 86 a does lead to diver- 
gence. This is surmounted by correcting the streamline position by an 
amount smaller than the calculated 6 n, or by utilizing an additional 
factor in Equation 86 a where is a constant which is greater than or 

equal to unity. 


A^(SA) D^6n) ^ - Aq ) 

AW 1 Ds^ AW 


( 86 b) 


The effect of setting to a value larger than one, say 1.5 or 2.0, is 
to reduce the curvature change between successive iterates and in an ap- 
proximate manner account for the items which were dropped from Equation 
82. 


5.6.2 Flow Inlet and Flow Exit Boundaries 

Options are provided for two types of boundary conditions at the upstream 
and downstream field boundaries. Either the flow angle can be specified or 
the curvature can be specified. The recommended boundary condition for 
general usage is that the curvature is zero (or, equivalently, the static 
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pressure is constant). In either case. Equation 86 a is applicable. However, 
for an upstream boundary, the coefficients Gi, G 2 and G 3 in Equation 85 
will all be zero. That is, the curvature correction at the central point 
will be related to the downstream points 4 and 5 only. If a curvature 
boundary condition is employed, then the curvature correction necessarily 
will be zero. Hence, all the G coefficients are zero, and the equation is: 

^2(6A) _ A2(Ax - Ao) 

AW AW 


The solution to this equation is simply: 

6 A = Ax - Ao (88) 

5,6.3 Body Surface Points 

The correction equation for the grid points on a body contour is 
trivial. The streamlines are already correctly positioned so: 

6 n = 0 (89) 

for each such point, 

5*8,4 The Far^-Field Interface Streamline Correction Equation 

In this section the correction equation for a constant pressure or 
far- field boundary is formulated. The notation is similar to that of Section 
5.6.1 and is illustrated in Figure 32, 

Again we begin with the continuity equation. The desired solution is 

(*H - Vl) = *N - »ll-l <»»> 

But for the velocity, V^, associated with the ‘Assumed** fai^field boundary. 


the equation is: 

(*XH - = «» - »N-1 

Equation 91 is subtracted from Equation 90, 

(*H - *H-iKpV>K.J - (A*)) - *xN-l) <Po’'o>B-J = » 

The identities: 

“ ^N-1 ^ (^oN " ^oN-l) “ ^^N-1 (93a) 

^xN "" “^xN-l “ (^oN ‘^oN-l) (“^x ” "^o) N " ^ ^o)n- 1 (93b) 
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Figure 32. Notation for Section 5.6,4. 


4 



and: 




(93c) 


are introduced into Equation 92 and the result rearranged to give: 
|<PoVc>».J + 6<pVVj) (6A* - 6A^.j) (A„„ - A„^J 5<pVYj 

= <PO^O>jj_J [(Ax “Aq[jJ -(Ajj - Aq)jJ_^[ 


Equation 94 represents the continuity form of the correction equation 
for the outermost streamtube. 

The velocity correction is approximated as follows. The average 
velocity, for the streamtube is related to the boundary velocity 

by the identity: 

^ + 6(Vj,_^ - Vn) (95) 

For a constant pressure boundary 6 Vj^ is zero; for the far-field interface 
streamline: 



5 

i=2 ®FFi ®"i,N 



(96) 


(97) 


Equations 96 and 97 relate the change in velocity at point (i = 4) to the 
local ordinate changes of the streamline (6nijN), Included in this for- 
mulation is the contribution from the movement of the point itself and the 
two upstream and two downstream neighbors. The coefficient (0.865/AS) and 
the relative factors of [0, -1, 2, -1, 0] were determined (as reasonable 
approximations) from numerical solutions. AS is the average local spacing 
between points along the interface streamline, and the coefficient (1-M§)“z 
corrects for the compressibility or Mach number effect. 

The second term on the right-hand side of Equation 95 is evaluated, 
as before, by the cross- stream momentum equation; 


1 

V 


9n 


= -c 
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Approximately} for the constant stream tube: 


^N- 


-i (n^ - tiN-l) (Cn + Cn-i) 


( 98 ) 


Define Vjj_^ = i (Vj^ + and Cu_^ = i (C^ + %-i>» 

can also be expressed as follows: 

(^N-i " ^n) “ ^ (“N ■ “N-i) S-i 


then Equation 98 

(99) 


The differential (or small variation form ) of Equation 99 is: 

■ %) = i (“n " °N-i) ^ ” 

= i - “N-l) 

^ (pN-i Vi (“n ■ “N-i)| ■ 


ViK vi ) 

PN-i 


^PN-i 


«N-l)l 

( 100 ) 


From the continuity equation, the first tern in the braces is zero. (For 
the outer streamtube, the effect of a change in radius is neglectable. ) 
The density variation in the second term can be related to the velocity 
variation: 


PN-i 



6 V, 


N-i 


Vi 


(101) 


Thus, Equation 100 can be written: 

^(Vi “ %) = i ^^(nN - tiN-i) 6Ci^i + iM^i CN-i(nN - Nn-i)6Vj^^ (102) 

Equations 96 and 102 are substituted into Equation 95 to obtain: 

^oN 5 

^-i ^Vi = 8^ i?2 Vi ^“i,N ^ Vi(“N " Vl)^Vi 

where: 

ON-i = 1 " i^N-i Vi (“N " ”n-i) 
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The product of curvature and streamtube thickness is expected to be very 
much less than unity, therefore Qu-J is approximately equal to unity. 
Also: 




(105) 


This together with Equation 103 substituted into Equation 94 gives the fol- 
lowing for the interface streamline correction equation: 


r y g 

(p^)n-J (^^N-1 “ ■ (^oN ■ 1^2 ®"i»N 

* * Vj ("» ■ "K-l)'®*-*] = [{*' ■ ■ (** ■ *°)>'] 


Equation 106 presumes: 

Also, to a reasonable approximation: 

®H-i ^ ®N 
PH-i ^ Pn 
Qu-i « % 

^oN 


Equation 106 is divided through by (pV)jj_i; and, with the above approxima- 
tions, the far- field boundary correction equation becomes: 

^N-1 ■ - <AoN - AoN- 1) — ®FFi ^"i.N " 

= <A* - Ao)h- 1 - (Ax - Ao)n (107) 

where: 

®N = ■ ^oM^l) ®N (“N - “N-l) <108 

In the present STC code, we further approximate Qn as unity. 
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The average streamline curvature is taken to be equal to the boundary 
curvature: 

-6Cn_^ = -SCn = ®3^“3,N + V®4,N + ®5^“5,n ®6^“e,N 

where the G*s are the influence coefficients relating streamline point 
movement to the negative of curvature. Also, it follows that: 

^%-l "" = ^^^N-1 ^^N-1 “ ^TT^N^n (110a) 


for axisyiranetric flow, and that: 

6A|^^1 - 6Ajj = - 6n^ (110b) 

for two-dimensional flow. 

With Equations 109 and 110, the left-hand side of Equation 107 involves 
only the streamline correction quantities, 6n, and the right-hand side is the 
error term computed in the flow balance section. Equation 107 is used with 
Equations 86, 87 and 89 to obtain the matrix equation for 6i^ for M = 1, 

2 .... where M is the total number of grid points in the field. 

5.6.5 The Curvature Influence Coefficients 

To complete the formulation of the system of correction equations, 
it is necessary to relate the change in curvature to the movement of the 
streamline. Specifically the values of G in Equation 85 (repeated below) 
are desired. 

-6C4^k = Gl6“l,k + G2fin2,k G3^“3,k + + Ggdng 

+ = 2:0^6 Hi 


An illustration of the notation is provided in Figure 33. To a very good 
approximation, the curvature variation is equal to the variation of the 
second streamwise derivative of 6n. 


-6C 


4,k 



(Repeat of Equation 86) 


To evaluate this second derivative, the linear spline equations are employed, 
maintaining compatibility with the beam curve fits of Section 5.2. The 
method for calculating the spline influence coefficients is presented in 
Appendix B. Since only small adjustments from the given curve are to be 
made, the arc length along the curve, s, and the normal point adjustment. 
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Assuned k 
Streamline 


4C, is the change in 

curvature between the 
assumed and corrected 
streamline at this point. 


Figure 33. Notation for Section 5.6,5 



6n, are replaced by x and y in that development. The resulting values of the 
influence coefficients are tabulated in Table II for the specific case of a 
uniform spacing where = As = 1. 

As indicated in Table II, the influence coefficients for a point on a 
spline curve involve all of the other points on the spline. That is, the 
movement of any one point affects the curvature at all of the points, al- 
though the effect on the curvature at distant points is small. To include 
the entire set of spline influence coefficients is Impractical, because it 
would lead to a completely dense, rather than sparce, coefficient matrix for 
the correction equation. Therefore, the curvature influence coefficients are 
always truncated to 5 points, or less near boundaries. The method of trun- 
cating is to simply replace the one long spline with a series of short 
splines. Each short spline passes through only the point in question and 
the two upstream and two downstream neighbors, if they exist. If the local 
velocity is supersonic, then three (or two) upstream points and no downstream 
points are included to duplicate, in effect, the supersonic curve fit of 
Section 5.2.3. 

Artificial boimdary conditions are enforced at the upstream and down- 
stream ends of this truncated spline. For subsonic flow, the selected 
boundary condition is that y" = 0 at each end (remember y = 6n). However, 
if the truncated spline end point is also a streamline end point, then the 
boundary condition is chosen to agree with the curvature end condition dis- 
cussed in Section 5.2.2. For supersonic points, the end options are again 
chosen to be equivalent to those used for the calculation of curvature. 

It is interesting to notice how different the spline coefficients are 
frcHn those of a pol 3 rnominal . This is illustrated in Table II. Although 
the 3-point parabola formula is commonly used for subsonic flow, the spline 
formula is found to be considerably different and, we believe, more 
accurate. 

5.7 MATRIX SOLUTION PROCEDURE 


In the previous section, equations for the streamline position cor- 
rection, 6n, were developed for each grid point. At Interior points and 
flow inlet/exit boundary points. Equation 86 is used together with 
Equations 84 and 85. And either Equation 89 or Equation 107 through 
110 are used at the orthogonal line boundary points. This yields a system 
of simultaneous equations where at any "central" point the applicable 
equation can be written in the form: 

Ai6i + ^ 3^3 ^4 J ^4 ^ ^ 5^5 ^ ^ 6^6 ^ ^ 7^7 ^ ^ 8^8 ~ 

( 111 ) 

The notation 6 is used in place of 6n for brevity and the subscripts refer 
to neighboring points as indicated in Figure 34. If the velocity at point 
4 is subsonic, then point 1 is not included (Aj^=0) ; if the velocity is 
supersonic then points 5 and 6 are both omitted (Ag = Ag = 0). 
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Table II. 


































Figure 34. Arrangement of Neighboring Points for the General 
Subsonic/Supersonic Star. 



It should be noted that points 1 to 6, as shown in Figure 34, never 
include any of the special points which terminate a partial orthogonal. 
These points are skipped by extending the star, as required, to the next 
point in either the downstream or upstream direction. 

In all cases, Aj^, A2, Ag, A^^\ Ag, and Ag are the curvature influence 
coefficients times the coefficient Bi. For example: 

Ai = B-j^ 

A2 = Bi Gg 


n r 

A4 = G4 


etc. 

( 2 ) 

The values of A^, Ag, and A^ , for interior points, are related to the 
flow difference between streamlines. These coefficients, given by Equation 
84 for axisymmetric geometry, are: 


A 

7 “ w. - w, 


k-l 


2n r, 


k+1 


^8 W, - W, 
k+1 k 


A 


( 2 ) 

4 


-2tt r, 
k 




( 112 ) 


For a two-dimensional configuration, the (2n r) factors are replaced 
by unity. Double streamlines (which separate two adjacent channels) are 
disconnected; the finite difference equation of the form of Equation 111 is 
only written for the second of the two coincident lines. For the first line, 
the following equation is employed: 

^4b “ ^4a “ ~ “o^4b “ (“x " ”o^4a 

The "4b" and "4a" subscripts indicate the points which belong to the 
lower and upper channel, respectively. For noninterior points, the values 
of Aj to Ag are similarly defined according to the equations of the previous 
section. 

The set of equations defined by Equation 111 are solved by a "block" 
relaxation method. In this procedure, an initial guess for the values of 
6 is successively relaxed until the solution does not change; viz, the 
solution is constant within a specified tolerance. The options available 
for sweeping the field are: 
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a) Solving orthogonal line blocks, sweeping upstream to downstream 


b) Solving streamline blocks, sweeping from the centerline to the 
outer radius (or from small y to large y) 

c) Alternate use of a) and b) 

Included In the relaxation procedure Is an acceleration factor, p, which 
serves the same purpose as the overrelaxatlon factor In the point- relaxation 
method. For the two types of block relaxation. Equation 111 is rewritten 
as follows: 

Orthogonal Blocks : 

*7*7 [*f ’ + f>*4“]*4 * *8*8 = - [aj «i 7 7 A 363 7 

<l-p) a‘“ 6 ^ 7 A 563 7 A 363 ] (113a) 


Streamline Blocks: 


^ 1^1 ^ 2^2 ^ 3^3 [^4 ] ^4 ^ 5^5 

( 2 ) 


6°6 


= RHS 


A 767 + ( 1 -p) A^ 64 + Ag 6 g 


(113b) 


The terms Involving the streamline adjustments on the modified right- 
hand side of each equation are evaluated using the results of the previous 
Iterate. Orthogonal line blocks require the solution of a tridiagonal 
matrix, while streamline blocks require the solution of a pentadlagonal 
matrix. In both instances, the solution Is obtained by a left-right de- 
composition of the coefficient matrix followed by back- substitution 
The solution is assumed to be converged when the maximum change from one 
iteration to the next is within a specified tolerance. Specifically it 
is required that: 

Max|A 6 | ^ TOLRL * Max | 6 | 

where A 6 is the change in 6 between two successive sweeps at a given point. 

The acceleration factor, p, may be either t^ken as constant or allowed 
to vary with the sweep number. Using the results of Peaceman and Rachford 
(ref. 14), the acceleration factor is determined as a function of the sweep 
number and the total number of points in the field. The latter parameter 
is tentative, since the STC program uses a nonuniform grid, and the 
Peaceman-Rachf ord relation was developed for a uniform grid. 
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(114) 


_ . 2 nn 

p = PB + 2pA sin 


where : 

n s sweep number 

Pg = base acceleration factor 

p^ = half of the amplitude of the sinusoidal variation 

mi = number of grid points in the field 

The standard mode of operation is to alternate between the orthogonal 
and streamline block sweeping. When the sweep number reaches the stage 
where n = 2 a,^W/tt then n, as used in Equation 114, is reset to 1. The re- 
sulting variation in convergence factor is depicted in Figure 35. 

As in the case of successive point overrelaxation (SOR) , an optimum 
value of p corresponding to the minimum problem convergence time, may be 
determined. Selection of the optimum acceleration factors is discussed 
below. 


An internal flow test case with a 0.09 was chosen to optimize the 
acceleration factor for the alternating direction solution procedure. The 
matrix solution time vs. p was established by systematically reducing pg 
fr<xn unity to a minimum value of 0.45. The results are depicted in Figure 
36. 


As indicated, the optimum p£ setting occurs somewhere in the vicinity 
of 0.5 to 0.6. With PB = 0.4, the block relaxation diverged. For pg = 0.5, 
the matrix solution time was 42.4 percent of the total iteration cycle time 
which includes one pass through Steps 3 through 12 of Section 3.0. 

Based on this test and other results similar to this, Pb and have 
both been initialized to 0.5 in the computer code. Alternate values may be 
input by the user. 

5.8 STREAMLINE ADJUSTMENT 

To adjust the streamline position, the coordinates are moved in the 
normal direction by the computed 6n's. 

z = z - Cy on sin 0^^ 

r^^ = r** + 6n cos 0^^ 
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Solutions , Seconds 


Sweep Number, n 


Figure 35. Typical Variation of the 
Acceleration Factor, 


GE635 Processor Time 
for 756 Grid Points 



Figure 36. Matrix Relaxation Time as a Function 
of the Base Acceleration Factor. 



The streamline angles, 0^^, are those found from the curve fits of Section 
5.2. The superscript p is the iteration sequence number, and is a conr 
vergence factor which is generally equal to unity. 

The new values of z and r then form the basis for the next iteration, or 
for the solution which is printed out. The next step in the procedure is to 
again compute the curvatures and from these compute the velocities and flow 
balance errors. 

5.9 BOW SHOCK WAVE 

The original Intent of the analysis development was to Include the 
approximate bow wave method of Moeckel (Reference 15.) With the inclusion 
of the bow wave, it would have been possible to analyze nacelles at transonic 
Mach numbers above 1.0. During the development of STC with the bow wave, 
serious difficulties were foimd and a reliable computer solution was not 
practical. Accordingly, two versions of the STC computer program were fur- 
nished; 1) a standard, fully functional STC computer analysis without the 
Moeckel bow wave and 2) a status level program incorporating the bow wave 
which was provided for the purpose of future investigation. 

Technically, the approximate nature of the Moeckel method compromises 
the quality of the STC solution. The assumed hyperbolic shock shape and the 
empirically specified stand-off distance force continuity to be violated 
between the shock and the body. In addition, once the Moeckel bow shock 
is added to the flow field (after three initial flow field refinements to 
stabilize the streamlines and orthogonals) , further refinement of the flow 
field in the vicinity of the shock is stopped. At the same time, streamline 
curvature near the shock is preset to zero. The good features of STC are 
not part of the bow wave solution. Hence, the STC results are not enhanced 
by including the approximate bow wave solution. 

The STC solution with the bow shock has been programmed as a separate 
computer program since extensive changes were necessary. The difficulties 
that developed during checkout preclude a useful solution, but Volume 11 
of the User's Manual (Reference 16) shows how to use this version of STC. 

For Mach numbers up to 1.2, it is recommended that the standard version 
of STC be run in the supersonic mode with grid refinement suppressed near 
the shock. At Mach 1.2, the nominal streamline deflection through the shock 
is 4° with a maximum total pressure loss coefficient of 0.007. These con- 
ditions can be approximated using the standard STC program. 

5.10 INTEGRAL MOMENTUM CHECKS AND PRESSURE DRAG EVALUATION 

The STC Program evaluates the thrust/drag on each boundary surface and 
then verifies these forces by performing "overall" momentum balance checks 
for each of the fluid streams. For a typical nacelle configuration, a 
momentum balance is computed for the inlet stream, the external stream, and 
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the jet exhaust stream. For the inlet stream, as an example, the entering 
momentum (or ram drag) calculated using fluid properties at the most upstream 
station is first determined. Second, the pressure times the projected area 
of the cowl stagnation streamline is computed. This is added to the axial 
pressure force on the underside of the inlet lip from the stagnation point 
to the last calculation station (near the fan face). Third, a similar 
pressure-area integration to include the axial forces on the spinner is 
performed. Finally, the sum of the entering momentum and the integrated 
upper and lower boundary forces (including the additive drag of the cowl 
approach streamline) is compared to the integrated axial momentum flux at 
the last station inside the inlet. A discrepancy will indicate inaccuracies 
in the computed pressure distributions or, perhaps, insufficient refinement 
of the calculation grid for adequate resolution. It has been found that 
these momentum checks are quite valuable for quickly assessing the computed 
result. 
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6.0 CONCLUSIONS 


The Streamtube Curvature Analysis has been developed and a computer 
solution has been utilized to solve the transonic flow field over an iso- 
lated nacelle. The analysis includes the capabilities to: 

• Analyze the stagnation region of an inlet with grid refinement 

as needed so that the stagnation streamline is properly predicted. 

• Handle multiple streams of differing stagnation properties with 
a static pressure balance at the interface. 

• Predict the location and strength of imbedded shock waves on the 
external nacelle surface. 

m Define an analytical far-field boundary so that free-f light con- 
ditions are predicted. 

• Achieve greater solution economy for the transonic flight speed 
regime in that the computational times are 5 to 10 times faster 
than state-of-the-art time- dependent methods for the same number 
of grid points. 

• Provide a user-oriented design analysis tool for responsive 
solutions to engineering problems. 

The inviscid pressure distributions predicted by the Streamtube Curva- 
ture Analysis compare very well with test results when viscous interactions 
are minimal. When shocks are predicted, the viscous interactions and bound- 
ary layer separation regions prevent good correlations with test data. 

The computer analysis fills a void in that several techniques are 
available for solving the inviscid equations of motion about arbitrary two- 
or three-dimensional bodies at transonic speeds, but none have handled the 
complete nacelle with inlet flow and exhaust flow. The Streamtube Curvature 
Analysis provides the capability to predict the transonic flow field about 
typical aircraft engine installations in isolated nacelles at transonic 
speeds. The solution technique provides a design analysis tool which will 
provide guidance for wind tunnel testing to develop nacelle shapes to min- 
imize drag within given design restraints. 

The details of the computer program operation, usage, and structure 
are documented in the Users Manual (Reference 16). 
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7.0 APPENDIX A 


SEC0ND--0RDER STREAMLINE CURVATURE EQUATIONS FOR 
ISENTROPIC PLANAR FLOW 


In this section the second-order equation which describes the shape 
of the streamline is derived. An orthogonal system is chosen. Hie variable 
n is a measure of the distance across the streamlines and s is the distance 
along a streamline. The partial differential operator 9 indicates that the 
direction of differentiation is normal to the streamline. The operator D 
indicates that the direction of differentiation is along a streamline. The 
basic equations are: 

Continuity: 

Crocco Form of Normal 
Momentum Equation for 
Irrotational Flow: 

where : 

V ss velocity 
p s= density 


C = streamline curvature =s - 

s s= streamwise coordinate 
n =s cross-stream coordinate 

0 = streamline angle measured from horizontal 

Equation 115 is differentiated as follows: 


pv an « aHf 


-cv 

9n 




(115) 

(116) 



(117) 


i /i i lfi\ ^ is 
~ pv \V p av/ 5n b'H 


(118) 
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Equation 117 is substituted into Equation 118 and the result rearranged: 




(pV)2 V p BV; V Bn 


(119) 


Now Equation 116 is substituted into Equation 119: 


BY^ 




(pV) 




( 120 ) 


Equation 120 is the desired second-order equation for n, where the curva- 
ture is analogous to a derivative of the form (D2n/Ds2) . Since the flow is 
isentropic, the coefficient in parenthesis may be expressed as follows: 


1 



1-m2 


( 121 ) 


thus, the equation reduces to: 
^ — 2 ' = » 


( 122 ) 


Clearly, the elliptic and hyperbolic nature of the equation is evident 
when the Mach number is less than and greater than unity. 

Equation 120 is now written in the following abbreviated form: 






2 


B C = 0 


where: 


B = 


1 

(pV)^ 




(123) 


The equation in its present form cannot be applied directly to obtain 
a solution for n because of the difficulty of relating the radius of curvature 
to the second streamwise derivative of n. Instead, it is used to calculate 
the streamline adjustments (in the cross-stream direction) for an assumed 
streamline pattern* 
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To illustrate, we consider an assumed set of streamlines which pass 
through the small circles in the adjoining sketch. 



By a curve-fitting process, the values of streamline angle and curva- 
ture, Cq, are determined, (The circle points are also orthogonal ized. 

That is, they are moved along the streamlines so that the ’’orthogonal 
lines” are truly normal to the given streamlines.) Equations 115 and 116 
are integrated along the orthogonal lines, assuming the value of is 
valid, and from this ’’flow balance” the x- positions of the streamline are 
determined. We have then the following equation satisfied: 



B Co = 0 


(124) 


The ”o” subscript denotes values related to the assumed streamline 
positions and curvatures. The ”x” subscript refers to the streamline posi 
tion as calculated by the continuity equation. The true solution to 
Equation 123 is sought; ’’true solution values” are unsubscripted. 

Equation 124 is subtracted from Equation 123, The result is: 



The adjustment to be made 
the above can be written: 

+ b(-C + C \ = 

\ ° / 


Finally, we note that: 


in streamline position is 6n = 


Hence 


(125) 
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(126) 



So the differential equation becomes: 


5^(6n) ^ „ D^(6n) 

Ds^ 


2 (®x~®o^ 


( 127 ) 
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8.0 APPENDIX B 


SECOND DERIVATIVE INFLUENCE COEFFICIENTS FOR A SPLINE 


For the interval between points i and i+1, a cubic equation may be 
written in the form: 

y = Yi + bi (x - x^) + 1/2 Ci (x - x^)2 + 1/6 (x - x^)? (128) 


Alternately, this same cubic may be expressed in terms of the ordinates 
(y^ and y±^i) and the second derivatives (c^ and at the points Xj^ and 

namely; 


y = 6^ [ci(xi+l-x)^ + Ci+i(x-Xi)3j 


(129) 





(x 


1+1 


X) 





(x-Xj^) 


where ~ *i‘ ®y requiring a match in the 'first derivative at the 

interval boundaries. Equation 129 may be used to derive a set of N-2 
equations for the c’s: 


AXi_i Cj^_j + 2(Axj^_l + Ax^)cj^ + Axi c^^j^ 


_ 6 i / 6 6 \ 6 

Axi_l yAxi_i * AXj^ * Ax^ yf+1 


(130) 


Two additional end conditions then yield a total of N equations. A 
variety of end conditions is possible. The following are available as 
standard options: 


(1) The first derivatives yj^' and/or yjj' are zero: 

6 0 

2AX1C1 + Ax^C 2 = - ^ yi + yg 


Axjj_iCN-i + 2Axn_iCn 


6 

Axjj_]^ 


6 

^N-1 - A^_^ yN 


(131) 


(132) 
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(2) The second derivatives yi" and/or yj^" 


are zero: 


(133) 


Cl = 0 


Cj, = 0 


(134) 


(3) The third derivative of the end polynomlnal is a specified 

fraction, F, of the third derivative of the polynomlnal second 
from the end: 


1 /I F \ F 

Lxi Ax2/ ^X2 “ 

^* N -2 (^* N -2 ^* N -1 ) ^^1 ^ 


(135) 

(136) 


(4) Same as Option 3 except that yj^ is not chosen arbitrarily. 
Instead it is required that yf^t = ~oyif (Note this option 
is coded only for the downstream end of the spline-) 


Axn-2 


* °^Vi ) 


®N-2 


4«h-2 I* * * 


^1 + oAxn-iJ 


CN-1 


AXjj_ 


N-2 


r 1 a ' 

^N-2 ® Axjj_2 * 1 + Oi^*N_i 


FN-1 


(137) 


Equations 130 through 136 can be written in matrix form: 

Aj.i ®i = ®j,i yi 


(138) 


where A and B are square coefficient matrices, and c and y are column 
vectors of length N. If Option (4) is used, then Equation 137 replaces 
Equation 130 for i = )f>l, and the size of the arrays are reduced by one. 
Generally, splines are fitted to given sets of data and Y is known. Here, 
we find the influence coefficients G by premultiplying Equation 138 by 
the Inverse of A. 
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cj = Gj,i yi 


(139) 


where : 

G = A"^B (140) 

G is the array of coefficients which give the influence of i = 1, 2..,N 
on the second derivative cj. 

We comment that G, as used in this section, is a square matrix. The 
desired set of influence coefficients for the second derivative at point j 
is the row of the matrix. In section 5.6.5, the notation Gj^ is used to 
indicate a vector which is the row of the matrix of this section. 
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9*0 APPENDIX C 


NOMENCLAIURE 


Symbols 

A 

A 

vv 


a 

B, Bj, 




b 

C 


c 



F 

f 

g 

G 

H 

M 

M 

o 

n 


flow area measure normal to the streamline = JaiTrdn 
coefficient matrix 
mass flow ratio 

unperturbed interface streamline radius 

coefficients in the streamline correction equation 

first derivative 

curvature, = - d0j^/ds 

second derivative 

inlet drag coefficient 

pressure coefficient 

specific heat at constant pressure 

integrated pressure drag coefficient 

diameter ratio 

streamline curvature curve-fitting parameter, see Section 5.2*2 

fractional position in the interval 

1-f 

curvature influence coefficients 

enthalpy 

Mach number 

free-streara Mach number 

distant measure along an orthogonal 
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n 

P 

Q 

q 

R 

g 

r 

S 

s 

T 

u 

V 

V 
W 

' Max. 

z 

a 

a 

3 

p 

V 

c 


n 

? 



matrix relaxation sweep number 
pressure 

parameter defined by Equation 104 

dynamic pressure 

perfect gas constant 

radius 

entropy 

distance measured along a streamline 
temperature 

axial component of velocity 
vertical component of velocity 
total component of velocity 
cumulative flow rate 
normalized inlet length 
rectangular coordinate axis 
axial position 

angle of the chord between two adjacent points on the 
curve 

parameter relating velocity corrections on the two 
sides of the jth slip line 

compressible similarity parameter, (1-M^)^ 

ratio of specific heats 

dummy z- variable 

dummy y-variable 

dummy x- variable 

coordinate system in the streamline direction 
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p 

p 

a 


0 

$ 


0 ) 


coordinate system in the direction normal to streamline 
fluid density 

matrix relaxation acceleration factor (Section 5«7) 

source density 

streamline angle 

orthogonal line angle 

velocity potential 

stream function 

vorticity 


Subscripts 

a 

a 

b 

b 

FF 

i 

j 

M 

0 

GO 

s 

S 

T 

X 

1 


first end of the interval 

evaluated above a slip line 

second end of the interval 

evaluated below a slip line 

far-field interface streamline 

orthogonal line index (Section 5*5 & 5*6) 

streamline index 

field point index 

initial or assumed values 

free stream 

static 

unperturbed interface streamline position 
total 

calculated by using the assumed streamline curvatures 
streamline direction 


2 


orthogonal line direction 


Superscripts 

(1) Denotes streamwise connection 

( 2 ) Denotes crosswise connection 

* 1st derivative 

»• 2nd derivative 

3rd derivative 

p iteration sequence number 
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10.0 ADDENDUM - PART II, TURBULENT BOUNDARY LAYER AND TURBULENT SEPARATION 
PREDICTION METHODS, By D.J. Lahti and P.H. Heck, General 
Electric, Aircraft Engine Group, Cincinnati, Ohio 45215 


10.1 INTRODUCTION 


The current state of the art of computing the development of turbulent 
boundary layers on arbitrary axisymraetric and/or planar bodies does not allow 
for the complete solution of the governing equations without some assumptions 
being made. These assumptions define the relationship of the fluctuating 
quantities to the mean flow quantities in determining the turbulent transport 
properties. As a result, many computational methods have evolved over the 
past few years, each method depending on the closure assumptions and, to some 
extent, the limited data available to support those assumptions. Without 
enough valid data to positively substantiate or disprove these closure as- 
sumptions it becomes difficult, if not impossible, to assess their true 
merits. 

Although there is general wide-scale agreement that the finite difference 
methods do possess several advantages in computing turbulent boundary layers, 
they also possess some disadvantages. Hie primary disadvantages from the 
practical or engineering point of view is their relatively long computational 
times and the attendant higher costs. 

There is general agreement within the General Electric Aircraft Engine 
Group that not all problems require the detailed solutions provided by the 
finite difference boundary layer methods. A large majority of the design de- 
cisions made can be accurately and confidently supported by the less expensive, 
yet accurate, integral boundary layer solutions. As a direct result of their 
rapid computational times and demonstrated accuracy, they are very attractive 
for coupling with inviscid flow analysis programs. This coupled inviscid/ 
viscous computer program provides a very valuable and efficient engineering 
tool. Hie integral boundary layer method selected for coupling with the 
Streamtube Curvature Program incorporates the method of Stratford and Beavers 
(Reference 1?) and is discussed in the following section. 
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10. li TECHNICAL DISCUSSION 


10.2.1 Turbulent Boundary Layer Method 


Reference l8 examines several well-known integral methods for calcu- 
lating the boundary layer characteristics and correlates the results as 
functions of free-stream Mach number and Reynolds number based upon an equiv- 
alent flat plate length (see References 18 through 25). This examination 
finds that all the methods concur in the conclusion that the momentum thick- 
ness may be expressed in the form: 


e = f (M) X 

where: 


X 




(141) 


(142) 


For the flow on a flat plate, the Mach number, and thus the function P, would 

be constant so that X s x and 9 = f (M) x R^~ . In effect, then, all the 

methods agree in saying that the momentum thickness in a flow with pressure 
gradient may be obtained from the expression for a flat plate, provided that 
actual distance x is replaced by an equivalent distance X, according to 
Equation 142. Consequently, the methods can only differ from one another in 
the expression for the flow on a flat plate and in the value of the function P. 
The only part of the result which can differ due to differences in each of the 
analytical treatments is, therefore, the function P. 

Five of the methods employ transformations, each being some form of 
Stewartson's transformation from compressible to incompressible flow. There- 
fore, differences in P are due to differences carried over from incompres- 
sible solutions. Stratford and Beavers then conclude that rather than trying 
to assign any priority due to the merit of the incompressible solutions, they 
rather choose to represent P by a reasonable average of all the solutions. 

This is the concept employed by Thwaites (Reference 26) for incompressible 
laminar flow, and it proves to be the most expedient method under the 
circumstances . 

In summary, the calculation procedure is as follows. 

The boundary layer momentum thickness is expressible in the form: 

9 = f (M) X (143) 

where X is the equivalent flat plate length. It is defined as the length 
over which a boundary layer growing on a flat plate at the given Mach number 
would acquire the same thickness as the real boundary layer at that given 
location. 


102 



Stratford and Beavers propose the following equations 
visional working formulas from which the following integral 
be calculated. 


For 


1 X 10® < Rx ^ 1 X 


i2/,rt\"0.7 -1/5 


0 = 0.036 (i + MVio) 


-1/5 


6 = 0.37 


6*= 0.046 (1 + 0.8 


For 


1 X lo' < RX < 1 X 10* 


9 = 0.022 (1 + m 2 / 10 )“®*'^ XRj ^" 1/6 


- 1/6 


6 a 0.23 XRj^ 

6 *= 0.028 (1 + 0.8 

1 * 

/ 


where X = — 


P dx 


for planar flow 


and X = 1 


/ 


Pr dx for axisymmetric flow 


P = Lh/(1 + M^/5)]^ 


and a = 


1.25 for 10 


1.20 for ^ 10' 


as the pro- 
parameters may 

(l44a) 
(l44b) 
( l44c) 

( l44d) 
(I44e) 
(l44f) 
(145a) 

(145b) 

(146) 

(147) 
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Having these expressions, the distributions of 9 (x) and 6*(x) can be cal- 
culated, given the boundary layer edge pressure distribution. Once the dis- 
tributions of 9 (x) and f>*(x) are known, the integral momentum equation 
can be solved numerically to determine the local skin friction coefficient. 


That is: 



2 dU 
U dx 


(20 + 


6 *) 


2 dr 20 dP 
r dx ^ P dx 


(148) 


where , 


1 0 for planar flow 

1 for axisymmetric 


(149) 


Equations l44 through 149 are those employed in the boundary layer 
solution. 


10.2.2 Example Cases - Boundary Layer 

The simplicity of the above method is quite obvious; and, as discussed 
in the Introduction section, it provides accurate and reliable estimates of 
the boundary layer effects so long as the assumptions of an adiabatic wall 
and turbulent flow throughout are not violated. In particular, it is ideal 
for calculating the integral displacement thickness distribution along a 
surface when an inviscid/viscous iterative calculation is being performed. 

Several sample cases which demonstrate the ability of this integral 
method to predict turbulent boundary layers are discussed in this section. 

In the examples which follow, the predictions of the Stratford and Beavers 
method are designated S-B. In order to provide some comparison of relative 
accuracy of this method with a finite difference solution, some of the exam- 
ples were also analyzed using the General Electric Aero Boundary Layer 
Program. This program, designed B.L. , is very similar to the program de- 
veloped at NASA-Langley by Beckwith and Bushnell (Reference 27). The results 
of the B.L. solutions on the various example cases are shown to provide a 
reference from which the quality of the S-B solutions can be judged. 

10.2.2.1 Incompressible Flat Plate 

The incompressible flat plate data of Wiegart (Reference 28 ) is shown 
in Figure 37 along with the S-B and B.L. predictions. The agreement is good 
over the full range of momentum thickness Reynold’s numbers. Reference 28 
suggests that the first data point is at about the minimum Reynolds number 
for turbulent flow. 

This is further verified by the B.L. solution which shows this to be 
at about the same location as the downstream end of the transition region. 
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Skin Friction Coefficient 



10 ® 10 ® 10 ' 


Reynolds Number 

Figure 37. Incompressible Flat 



10* 2. 2. 2 Flat Plate with Pressure Gradient 


Reference 29 reports the results of a study to experimentally de- 
termine the Reynolds Analogy factor on a flat plate mounted in the diverging 
section of a supersonic nozzle. Although most of the data taken in this 
study contained the effects of wall heat transfer, measurements of the mo- 
mentum and displacement thickness distributions were made for the adiabatic 
wall case. These data are shown on Figures 38 and 39 along with the S-B 
predictions. The e3q>erimental pressure distribution started about 25% of 
the length of the plate from the leading edge. Therefore, the experimental 
pressure distribution was extrapolated to the leading edge, and the boundary 
layer was assumed to start from there. The predictions indicate good agree- 
ment with the data. 

10.2.2*3 Waisted Body of Revolution 

Figure 4o is a sketch of the waisted body tested by Winter, Rotta, 
and Smith (Reference 30). This work is usually recognized as one of the 
best available sources of compressible turbulent boundary layer data on a 
surface other than a flat plate. In addition, these data are frequently 
used as a basis for comparison for the various boundary layer prediction 
methods. Therefore, it serves as a standard reference from which the quality 
of these methods may be compared. 

Figures 4l through 43 show the S-B and B.L. predictions of momentum 
thickness, displacement thickness, and skin friction coefficient, respectively, 
for a free-stream Mach number of 0.597* In general, the agreement with the 
data is very good for the momentum thickness and displacement thickness dis- 
tributions of both S-B and B.L. However, at this Mach number, the skin fric- 
tion predictions tend to be somewhat higher than the data. The B.L. predic- 
tions are closer to the data between X/L =: 0.3 ^nd X/L = 0.7 » ^nd the S-B 
predictions are closer between X/L = 0.7 and X/L = 1.0. 

Similar predictions are made for free-stream Mach numbers of 1.4o4, 

1.7, and 2.0. These predictions are compared with the data in B'igures 44 
through 52. In general, it can be said that in all cases the agreement is 
very good for both the S-B and the B.L. predictions. For engineering pur- 
poses, the primary advantage to using the S-B method is its extreme simplicity 
and almost negligible expense as compared to the finite difference method. 

10.2.2.4 Supersonic Ramp - Adverse Pressure Gradient 

Reference 31 reports the results of measurements of the turbulent 
boundary layer on each of three supersonic compression ramps. Although not 
a large quantity of data was taken, the skin friction coefficients were 
measured. Predictions using S-B are shown in Figure 53 along with the test 
data. The agreement is seen to be very good in all three cases. 


106 



012 


.C0888888888 


(ui3)0 ~ ssau^im, iun:).uaiuo|j 


107 


Figu: 
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Figure 39. Flat Plate with Pressure Gradient, Displacement Thickness Vs. Distance from Leading Edge. 






Axial Distance — X/L 


= .597 

^ 6 
Re/^ = 6.46 X 10 /m. 

Figure 41. Waisted Body Revolution, Momentum Thickness Vs. Axial 
Distance. M = 0.597. 



Displacement Thickness-^ i*(cm.) 



M = .597 

® A 

Re/j^ = 6.46 X 10 /m. 

Figure 42. Walsted Body of Revolution, Displacement Thickness Vs. Axial 
Distance, M = 0.597. 
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Re/j^ = 6.46 X 10 /tn. 


Figure 43. Waisted Body of Revolution, Skin Friction Coefficient 
Vs. Axial Distance, M - 0.597. 
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Figure 44. Waisted Body of Revolution, Momentum Thickness Vs. Axial 

Distance, M = 1.404. 
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Momentum Thickness 9 (in.) 


Displacement Thickness — i*(cm.) 
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Axial Distance ~ X/L 

M = 1.404 
o _ 

Re/j^ = 1.31 X lOVm. 


Figure 45. Waisted Body of Revolution, Displacement Thickness Vs. 

Axial Distribution, M = 1,404. 
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Re/j^ = 1.31 X lo'/m. 


Figure 46. Waisted Body of Revolution, Skin Friction Coefficient 

Vs. Axial Distance, M = 1.404. 
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Axial Distance X/L 


M = 1.7 

o ^ 

Re/^ = 6.61 X 10 /m. 

Figure 47, Waisted Body of Revolution, Momentum Thickness Vs. Axial 
Distance, M.. = 1.7. 
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Figure 48. Waisted Body of Revolution, Displacement Thickness Vs. Axial 
Distance, M = 1.7. 
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Figure 50, Waisted Body of Revolution, Momentum Thickness Vs. Axial 
Distance. = 2.0. 
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Figure 52. Waisted Body of Revolution, Skin Friction Coefficient 

Vs. Axial Distance, M = 2.0. 
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Figure 53. Skin Friction Coefficient on 
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10.2«3 Separation Prediction Method 


In general, the prediction of turbulent boundary layer separation is 
best accomplished by performing a finite difference boundary layer solution 
over the given body surface. Separation is predicted at the point where the 
local skin friction coefficient becomes zero. Although this method provides 
reasonably accurate results, it requires a complete boundary layer solution 
which may not always be desired or necessary. 

In the past the only other alternatives to a complete finite difference 
boundary layer solution were; l) a complete integral boundary layer solution, 
or 2) the determination of some pressure gradient parameter which takes on a 
certain value at or near the separation point. The integral boundary layer 
solutions were known to give relatively large errors in skin friction co- 
efficient near the separation point; and, as a result, predictions of sepa- 
ration were unreliable and inconsistent. Subsequent improvements in the 
integral methods allowed more accurate boundary layer skin friction calcu- 
lations and^ thus, separation predictions; however, they were accompanied 
by increased cost and complexity. Hie simple pressure gradient parameter 
methods of predicting separation have been used with some success for some 
flows. However, they have not received widespread use because they also have 
been unreliable and inconsistent. 

Reference J2 reviews several methods for calculating incompressible 
turbulent boundary layer separation, and concludes that the Stratford method 
(Reference 33) quite satisfactorily predicts it. However, as is often the 
case, there is very little information regarding its usefulness for compres- 
sible flows. Therefore, the Stratford method was modified for compressibility 
and exercised on several compressible flow configurations where boundary layer 
separation was measured. In addition, General Electric* s finite difference 
AERO Boundary Layer program was also run on several of these configurations 
to establish some sort of reference from which the quality of the Stratford 
predictions could be judged. Finally, this method has been incorporated into 
the NASA-Langley version of the STC/viscous analysis program. 

Reference 33 presents the complete development of the basic theory for 
predicting the separation of an incompressible turbulent boundary layer. 

Only the basic ideas will be discussed here. 

The method postulates that the turbulent layer can be divided into two 
distinct regions. In the outer region, the pressure rise just causes a 
lowering of the dynamic pressure profile, but does not change its shape. 

The losses due to shear stresses within this region are assumed to be almost 
the same as those on a flat plate under the influence of the same pressure 
rise. In the inner region, however, the inertia forces are too small to 
overcome the pressure gradient, and the velocity profile is distorted. In 
this region the pressure forces are balanced primarily by the tranverse 
gradient of shear. 
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Basically, Stratford pieces these two regions together, describing 
the outer region using the similarity power-law profile: 




(150) 


and the inner region using the reduced boundary layer momentum equation: 


St dp 
dx 


(151) 


Utilizing Equation 150 and the definition of the stream function 
Y = /iidy, one can arrive at the following expressions for the outer 
region: 


Sy 

u = U 

o 


(^) 


and, 




n U 6 
o 

n + 1 



+ 1 /n 


(152) 

(153) 


(154) 


Integrating Equation 15I and utilizing Prandtl*s mixing length expression, 
one can arrive at the following expressions for the inner region when the 
shear stress is equal to zero: 


du j 

f 1 

dP 

1 

1/3 

dy ■ 1 

ipv 

dx 

Py 

) 

u=| 

4 

dP 


1/2 

Uv 

dx 

0 / 



y I - iE \ ^ 

" 3 Ippv dx ; 


(155) 

(156) 

(157) 
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The K is the Karman constant, and P is an empirical constant intro- 
duced to account for the fact that Equation 151 is only approximate near 
the wall and exact only at the wall. In addition, the mixing length ex- 
pression changes with increasing distance from the wall. Therefore, it is 
expedient to account for these effects by incorporating the single anpirical 
factor P, whose value is determined by experiment. 

At the interface between the inner and outer layers, u, du/dy, and ^ 
must be continuous. Therefore, equating them and performing some algebraic 
manipulations one finally arrives at the expression; 


(2n-4)/n 

ii) 



where C 

P 


If it is 


P-^Po 

1/2 PU^ 
o 

assumed that: 



1 - — f- for incompressible flow, 

o 


6 


X 


0.37 Re 


-1/5 

X 


P = 0.66 

K a 0.4l and n a 6 

then Equation I 58 can be written as; 

/ ^ \ V2 / (■ v-0.10 

= F(x) S' 0.35 


(158) 


(159) 


(160) 


This equation assumes turbulent flow throughout and is only appli- 
cable in an adverse pressure gradient. Vhen there is a region of favorable 
pressure gradient or the turbulent boundary layer starts at other than a 
leading edge, the surface length coordinate X is replaced by S = X + X . 

The distance X represents the virtual origin of the turbulent boundary 
layer at the minimum pressure point prior to the start of the pressure rise 
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or adverse pressure gradient. The virtual origin X is calculated by using 
Equation l45a or l45b as shown in Figure 54. Hence, the upstream history 
of the boundary layer is included in the definition of the separation 
parameter. 


According to the analysis, Equation l 6 o will predict separation when 
F(x) = 0.35 for incompressible flow. However, for a typical turbulent 
boundary layer flow, F(x) increases prior to separation and decreases after 
separation. Therefore, after applying this method to several flows with 
turbulent separation, Stratford observed that if the maximum value of F(x), 

a. is greater than 0.4, separation is predicted when F(x) = 0.4, 

b. lies between 0.35 and 0.4, separation occurs at the maximum value, 

c. is less than 0 . 35 1 separation does not occur. 

This method works quite well for incompressible flows as shown in 
Reference 32. However, in order that the method be useful for solving prob- 
lems of more current interest, it must be extended to account for compres- 
sibility. The approach taken here is one dictated by expediency and physical 
reasoning rather than mathematical rigor. Basically the philosophy of the 
approach is as follows. 


The physical ideas postulated for the incompressible case are expected 
to apply in the compressible case. Therefore, rather than quantitatively 
alter the formulation of the problem to account for variable density, ass\ime 
that Equation I 60 is still valid but with the following exceptions. First: 


C 

P 


1 


2 

ue 



o 


is replaced by 



and second, the 


critical range of the function F(x) is now different from that for incom- 
pressible flows. 

Now, the only remaining task is to determine the new critical range 
for F(x). This is accomplished in the following way. By computing F(x) 
using Equation 160 and the pressure distribution for a compressible separated 
flow, determine the value of F(x) corresponding to the measured separation 
point. Once this is done for a number of cases, one can determine the upper 
limit for F(x). The lower limit is determined in the same manner, however, 
this time the calculations are done for compressible flows with adverse 
pressure gradient where it is known that the boundary layer does not separate. 
The greatest difficulty encountered in trying to determine the appropriate 
critical range for F(x) is the general lack of good tcompressible separated 
flow experimental data. Even though this difficulty exists, it is felt that 
there is enough good data available with which to determine this critical 
range. Figure 55 was derived from the data of References 32 through 38 . 
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Separation 



Figure 55, Separation Functions Vs. Peak Mach Number 





It shows the critical range of F(x) for determining if separation occurs. 

The only region for which no data were found is that region between M as O 
and M = 0.35^« This region is bracketed by the dashed lines; however, 
it is^expected that, for peak Mach numbers within this range, reliable engi- 
neering estimates can be made using the assumed linear variations shown. 


10.2.4 Example Cases - Separation 

The Stratford method has been used to evaluate several e3q>erimentally 
derived pressure distributions to determine how well it predicted the mea- 
sured separations. Cases have been found between incompressible flows and 
compressible flows at Mach numbers as high as 4.92. Each of these cases 
will be discussed. 


10.2.4.1 Airfoil - Incompressible Flow 

The well-known airfoil-type body tested by Schubaur and Klebanoff 
(Reference 34) was selected to check out the program. First the S-B boundary 
layer program was run using the test pressure distribution. At the minimum 
pressure point, the equivalent flat plate length is obtained directly from 
the boundary layer output. Hiis then serves as the length to the boundary 
layer virtual origin which is required in the separation calculation. 

Figure 56 shows the airfoil, measured pressure distribution, and the pre- 
dicted and measured separation locations. The agreement is seen to be very 
good. The AERO Boundary Layer Program did not predict separation to occur 
on the airfoil. 


10.2.4.2 Forward Facing Step - Subsonic Flow 

The subsonic compressible data of Chapman, Kuehn, and Larson (Ref- 
erence 35) is ideal for determining how well the Stratford method is able 
to predict separation because the exact boundary layer origin is known. The 
measured pressure distribution and predicted and measured separation locations 
are shown in Figure 57* The agreement is seen to be excellent. 

10.2.4.3 Circular Arc Airfoil - Subsonic Flow 

The circular arc airfoil data of Reference 36 were analyzed in a manner 
similar to that of the Shubaur and Klebanoff airfoil. That is, the measured 
pressure distribution was used to compute the boundary layer origin at the 
minimum pressure point. This equivalent flat plate length is then input into 
the STRTFD (S-B) program . This procedure was used for each of the four free- 
stream Mach numbers shown in Figure 58. The fact that the predicted separa- 
tion locations are closer to the measured ones for the 0.354 and 0.663 free- 
stream Mach numbers is probably due to the some^at uncertain initial boundary 
layer conditions. Since the airfoil is just a bump on the wind tunnel wall, 
with suction and blowing slots upstream of the leading edge, the actual bound- 
ary layer height is not known as accurately as is desired. This affects the 
equivalent flat plate length which is calculated at the minimum pressure 
point on the airfoil, which, in turn, affects the value of the predicted 
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separation point. Even though these differences occur, they introduced 
only small inaccuracies. 

10.2.4.4 Flat Plate Shock Boundary Layer Interaction - Supersonic Flow 

Seddon (Reference 37) has studied the shock boundary layer inter- 
action on a flat plate. The mechanism for introducing the pressure gradient 
on the flat plate is the normal/oblique shock created by a second plate lo- 
cated above the first. The incident shock wave interacts with the plate 
boundary layer and causes a steep, but not discontinuous pressure rise. The 
measured pressure distribution and measured separation point are shown on 
Figure 59 along with the predicted separation point. Again the agreement 
is seen to be excellent. 

10.2.4.5 Wedge - Supersonic Flow 

An experimental study of the conditions necessary to promote boundary 
layer separation in the compression corner created by the intersection of 
a wedge and a plane wall at very high Reynolds numbers is reported in Ref- 
erence 38 . In this study the wedge/wall intersection was formed by a for- 
ward hinged plate in the wind tunnel floor. Changing the wedge angle was 
accomplished by simply swinging the hinged plate through the desired angle. 
There were many conditions simulated in this test, and only three were ana- 
lyzed with the Stratford method. Hiese were selected randomly from all the 
available cases. Figures 60 and 6 l show the measured pressure distributions 
along with the predicted and measured separation points for a free-stream 
Mach number of 2.95- The only difference between the two figvires is the 
wedge angle. Figure 62 shows these same variables at a free-stream Mach 
number of 4.92. Ihese cases were selected to demonstrate the ability of 
the method to predict turbulent separations even for very high Mach numbers 
and Reynolds numbers. One should also note in Figure 62 that the finite 
difference Boundary Layer Program tends to predict separation sooner than 
does the Stratford method. 

10 . 2.5 Data Comparisons - Axisymmetric Inlets 

The predicted pressure distributions on the NASA ATT inlet No. 8 
(NASA I- 85 -IOO) have been compared with the measured surface pressures re- 
corded on the through-flow nacelle during testing in the l6-foot tunnel at 
NASA-Langley. In Section 4.3» at the beginning of this book, the emphasis 
was on the predicted results from the Streamtube Curvature analysis. It 
was evident then that the boundary layer displacement effects and separation 
had a significant influence pear the leading edge of the inlet lip. Now, 
the viscous characteristics will be discussed in more detail. Two figures 
from Section 4.3 at the beginning of this book will be repeated with ad- 
ditional information noted. 

At a free-stream Mach number of M =0.8 and a mass flow ratio of 
0 . 81 , the streamtube curvature analysis with viscous effects included pre- 
dicts the surface pressures shown in Figure 63 . On the internal surface 
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of the cowl Lip, the comparison with test data is excellent and the effects 
of the boundary layer displacement thickness are not noticeable (see Figure 
15 ), On the external surface, the initial acceleration around the cowl lip 
is followed by a sharp pressure rise, and incipient separation is predicted. 
Once the separation point is predicted, the boundary layer displacement 
thickness becomes a constant from there on (the viscous analysis is not ap- 
plicable to flow in a separation region or reattachment). Hie streamtube 
curvature analysis oscillates in the region where the test data indicate 
a separation bubble has been formed. When the flow reattaches the predicted 
pressures and the measured pressures agree very well. Overall, the measured 
pressure force on the cowl outer surface was s - 0.039 and the predicted 

pressure force was ^ P 

Dp = -0.047. 


Ihe predicted and measured cowl pressures are shown in Figure 64 , at 
M =0.85 and a mass ratio of O.8I. The trends are very similar to those 
a^ M s 0.80 in that a separation bubble is still present but it is extended. 
Hie pressure comparison is very good on the surfaces where attached flow 
exists. 


The NASA No. 8 inlet was analyzed at M = 0.9 and a mass flow ratio 
of 0.883 with the viscous effects included. The flow ren^pined attached on 
the outer cowl surface. On the internal surface, the flow became supersonic 
locally and then shocked back to subsonic. Again, the comparison between 
measured and predicted surface pressures is excellent (Figure 65). Hie over- 
all integrated pressure force on the outer cowl was measured as -0.026 and 
predicted as - 0 . 031 . The streamtube curvature analysis, with viscous effects 
included, matches the test results more closely than the inviscid analysis 
alone. 


10.3 CONCLUSIONS 

The Stratford and Beavers integral boundary layer analysis has been 
integrated jwith the Streamtube Curvature inviscid flow analysis. The analysis 
has been linked so that the inviscid flow field provides the free-stream bound- 
ary conditions. Including pressure gradients, for the boundary layer analysis. 
The viscous analysis, in turn, modifies the nacelle geometry for boundary 
layer displacement thickness effects. The comparisons with test data and 
another finite difference analysis show that the method is an accurate and 
reliable predictor of turbulent boundary layers. 

The Stratford separation method, modified for compressibility, has 
been developed to predict incipient separation in the adverse pressure gradi- 
ent associated with the flow about a nacelle and that associated with embedded 
shocks. Its validity has been demonstrated by theory to data comparisons for 
two Reynolds numbers. This combined set of flow analysis methods provides a 
very useful and reliable engineering analysis tool. 
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Figure 64. Inlet Pressure Comparisons @ Mq = 0.8497 
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Figure 65. Inlet Pressure Comparisons @ Mq = 0.9007 




10.4 NOMENCLATURE 


Symbols 
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U 

U 
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Description 
Empirical constant 
Skin friction coefficient 



o 


Boundary layer thickness 
Boundary layer displacement thickness 
Karman constant; K = 0.4l 
Mach number 

Exponent in power law profile 

Pressure; also a function of Mach 
number (see Equation l46) 

Stream function 

Radius from centerline 

Reynolds number based on reference 
velocity and local x 

Reynolds number based on local velocity 
and local x 

Density 

Shear stress 

Boundary layer momentum thickness 
Velocity 

Velocity at outer edge of boundary layer 
Reference velocity 
Coordinate along wall 
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X 


Equivalent flat plate length 


y 


Coordinate normal to wall 


SUBSCRIPTS 


e Conditions at outer edge 

o Denotes free-stream or reference 

conditions 


T 


Stagnation condition 
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